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OD REDAKCJI 


Czytelnikom pierwszego zeszytu Nauk Matematycznych i Technicznych 
„Szczecińskich Roczników Naukowych” należy się krótkie objaśnienie, dla- 
czego są one bardzo potrzebne, a zarazem wielce opóźnione 1 niekompletne. 

Pilność potrzeby najlepiej uzasadniają fakty: Szczecin jako największe 
miasto Pomorza Zachodniego nie miał dotychczas samodzielnej oficyny 
wydawniczej ani własnego międzydyscyplinarnego ciągłego wydawnictwa na- 
ukowego. A przecież od dawna istnieje tu pięć wyższych uczelni akademic- 
kich i szereg prężnych pozauczelnianych pracowni naukowych, pracuje blisko 
350 profesorów 1 docentów oraz uczy się około 14 tys. studentów. W różnych 
krajowych instytucjach wydawniczych opublikowano stąd ponad 250 nauko- 
wych druków zwartych, o łącznej objętości około 1900 arkuszy wydawni- 
czych, nie licząc tysięcy krótszych publikacji naukowych w czasopismach 
specjalistycznych ani licznych tytułów dydaktycznych, informacyjno-sprawoz- 
dawczych i upowszechniających. W tych warunkach najtrudniej jest tu 
szybko publikować oryginalne prace badawcze średniej objętości, jakimi 
mogą być zwięźle zredagowane rozprawy doktorskie i habilitacyjne, albo też 
podstawowe szkice większych problemów, które dopiero po dłuższym czasie 
dopracowywania mogą stać się pełnymi monografiami. 

Dlatego też „Szczecińskie Roczniki Naukowe” chcą przyjmować do dru- 
ku tylko wartościowe oryginalne prace badawcze, streszczone do objętości 1 
ark. wyd., z dyscyplin reprezentowanych przez poszczególne wydziały Szcze- 
cińskiego Towarzystwa Naukowego: Nauki Społeczne, Nauki Przyrodnicze 1 
Rolnicze, Nauki Medyczne, Nauki Matematyczne i Techniczne oraz Nauki 
Morskie. W zasadzie każda praca powinna być najpierw przedstawiona na 
posiedzeniu naukowym odpowiedniego wydziału STN albo Komitetu bądź 
Komisji PAN, a przed przyjęciem do druku musi uzyskać pozytywne opinie 
dwóch recenzentów. 

Opóźnienie startu można ocenić na 25 lat i tłumaczyć tym, że od zarania 
odbudowy polskiej państwowości na Pomorzu Zachodnim tworzenie, organi- 
zacja i rozwój szkolnictwa wyższego 1 nauki dokonywane były na zgliszczach 
wojennych, w pustyni kulturowej, zupełnie od nowa, nieomal z niczego — 
głównie z heroicznego wysiłku ówczesnych pionierów. Działo się to bardzo 
szybko i wydajnie, ale chyba też nie bez pewnych omyłek. Szybki rozwój 
wyższych uczelni, głównie o dużym znaczeniu praktycznym, bez równoczes- 


nego powołania uniwersytetu, pozostał bez instytucjonalnego czynnika scala- 
jącego. W tych warunkach rolę koordynatora społecznego pełnić zaczęło 
Szczecińskie Towarzystwo Naukowe, które położyło w tym okresie wielkie 
zasługi, ale nie miało środków finansowych i mocy urzędowej, aby na 
przełomie lat pięćdziesiątych 1 sześćdziesiątych zorganizować bardzo już 
potrzebne własne ciągłe wydawnictwo naukowe. Możliwość taka zarysowała 
się dopiero w ostatnich latach, kiedy to dawne i usilne starania STN poparte 
zostały przez Wojewodę Szczecińskiego i Szczecińską Filię Oddziału Polskiej 
Akademii Nauk w Gdańsku. 

Podobnie też tylko z powodu ograniczonych środków finansowych drugi 
tom naszego wydawnictwa ukazuje się jako wysoce niekompletny — w 
postaci zaledwie dwóch zeszytów 6-arkuszowych, chociaż jest zaplanowane i 
powinno być drukowane co roku w postaci 5 zeszytów po 10 arkuszy w 
każdej sekcji naukowej, plus 5 arkuszy na zeszyt 6 poświęcony przeglądowi 
dorobku i kronice szczecińskiego środowiska naukowego za rok ubiegły. 

W końcu Redakcja pragnie wyrazić należne podziękowania: Władzom 
Województwa Szczecińskiego za zrozumienie potrzeb w zawsze trudnej i 
kosztownej działalności na polu nauki i kultury oraz ostatnie poparcie tego 
liczącą się już subwencją na cele wydawnicze; Dyrekcji Zakładu Narodowego 
im. Ossolińskich — Wydawnictwo we Wrocławiu za łaskawe przyjęcie Ro- 
czników pod zaszczytną dla nas opiekę wydawniczą oraz licznym Autorom 
za oferowane prace i Recenzentom za terminowe opracowanie ocen, zwła- 
szcza najtrudniejszych opinii negatywnych; do Czytelników zaś zwracamy się 
z uprzejmą prośbą o łaskawe nadsyłanie prac do druku, jak również kryty- 
cznych uwag o tomie wydanym i konstruktywnych propozycji na przyszłość. 


Eugeniusz Miętkiewski 
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ROMAN SOBAŃSKI 


TERMICZNE I KALORYCZNE EFEKTY MIESZANIA 
GAZÓW RZECZYWISTYCH 


Katedra Techniki Cieplnej Politechniki Szczecińskiej, 
Al. Piastów 19, 70-310 Szczecin 


Streszczenie. W pracy omówiono niektóre zagadnienia związane z miesza- 
niem się gazów rzeczywistych. Wykazano, że przy spotykanych w praktyce para- 
metrach (ciśnienia < 20 MPa i umiarkowane temperatury), mieszaniny wielu 
gazów należy traktować jak roztwory niedoskonałe. Procesowi mieszania towarzy- 
szą wówczas wyraźne efekty cieplne oraz zmiany parametrów termicznych. Zjawi- 
ska te występują między innymi podczas sporządzania mieszanin oddechowych 
stosowanych przy nurkowaniu saturowanym. Wyprowadzono wzory umożliwiają- 
ce obliczenie termicznych i kalorycznych efektów towarzyszących procesowi mie- 
szania gazów rzeczywistych. Przytoczone w pracy wyniki obliczeń dla mieszanin 
helowo-tlenowych wskazują na konieczność uwzględniania tych efektów podczas 
projektowania systemów urządzeń do badań podwodnych. 

Słowa kluczowe: termodynamika — mieszaniny gazów rzeczywistych. 


WSTĘP 


Przy niskich i umiarkowanych ciśnieniach mieszaniny gazów rzeczywi- 
stych zachowują się jak mieszaniny doskonałe spełniające prawa Daltona i 
Amagata. Wtedy 


p(V, T, Xx;) =2p,(V, T) (1) 
oraz 
(p, T, x,) = 2%; 5,(p, T). (2) 


W konsekwencji również energię wewnętrzną, entalpię oraz ciepło właściwe 
takiej mieszaniny można obliczyć na podstawie zasady addytywności ze 
wzorów 


a(V, T, x;) 22 u,(V, T); (3) 


h(p, dz; Xx,) => x;h;(p. T); (4) 


10 
c„(V, T, Xj) ej >, X;E(V, +); (5) 
J 


6, (DY TX;)= 2%; C„(P, T). (6) 


W miarę wzrostu cisnienia (zwłaszcza przy niskich temperaturach) coraz 
silniejsze stają się oddziaływania międzycząsteczkowe, w wyniku czego drobi- 
ny danego składnika w mieszaninie zachowują się inaczej (przy tych samych 
parametrach termicznych) niż w stanie czystym. Jedną z konsekwencji tej 
różnicy we wzajemnym oddziaływaniu pomiędzy cząsteczkami poszczegól- 
nych gazów w stanie czystym i w mieszaninie jest odstępstwo od praw 
Daltona i Amagata, a tym samym od zasady addytywności. Znajomość tych 
odchyleń, charakterystycznych dla mieszanin niedoskonałych, jest niezbędna 
przy rozpatrywaniu wielu procesów technicznych. Jednym z nich jest proces 
przygotowania i użytkowania hiperbarycznych atmosfer oddechowych stoso- 
wanych podczas nurkowania saturowanego, umożliwiającego człowiekowi 
długotrwały pobyt pod wodą. Zgodnie z fizjologicznymi wymogami organiz- 
mu, podczas nurkowania saturowanego człowiek oddycha mieszaniną gazów 
o składzie i ciśnieniu odpowiadającym głębokości nurkowania [3, 5, 7]. Poza 
tlenem, głównymi składnikami takich mieszanin są lekkie, obojętne dla 
organizmu gazy (np. hel, wodór lub neon). 

Gotowe mieszaniny lub ich składniki sprężone do ciśnienia przekraczają- 
cego często 20 MPa, magazynuje się w ciśnieniowych zbiornikach, a następ- 
nie wykorzystuje do sporządzania hiperbarycznych atmosfer oddechowych o 
dokładnie określonych parametrach. Od dokładności ustalenia parametrów 
przygotowanej mieszaniny zależy często zdrowie, a nawet życie oddychające- 
go nią nurka. Dlatego przy projektowaniu i eksploatacji systemów do 
głębokiego nurkowania niezbędna jest dobra znajomość wszystkich właści- 
wości stosowanych mieszanin, jak również termicznych i kalorycznych efek- 
tów towarzyszących procesowi mieszania gazów rzeczywistych. Efekty te 
można wyznaczyć eksperymentalnie lub analitycznie. Poniżej przedstawiony 
będzie analityczny sposób określania termicznych i kalorycznych efektów 
występujących podczas sporządzania hiperbarycznych mieszanin oddecho- 
wych. 


MIESZANIE GAZÓW RZECZYWISTYCH W WARUNKACH 
IZOBARYCZNYCH (p = idem) 


Proces izobarycznego mieszania gazów można prześledzić na modelu 
przedstawionym na ryc. 1. W izolowanym termicznie pionowym cylindrze ze 
szczelnym tłokiem, który pozwala na utrzymanie stałego ciśnienia, znajdują 
się gazy służące do sporządzenia mieszaniny gazowej. Są one przedzielone 
nieprzepuszczalnymi, diatermicznymi przegrodami. W stanie równowagi ter- 
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Ryc. 1 


modynamicznej ciśnienia i temperatury wszystkich gazów są jednakowe 
(ryc. la). 

Przy usunięciu przegród gazy mieszają się adiatermicznie ze sobą, two- 
_ rząc mieszaninę o parametrach jak na rycinie Ib. Oddziaływanie pomiędzy 
cząsteczkami składników mieszaniny prowadzi zwykle do zmiany średniej 
odległości pomiędzy drobinami, a tym samym do zmiany objętości gazu w 
cylindrze. Przy wzroście objętości rośnie energia potencjalna drobin, co w 
warunkach adiatermicznych pociąga za sobą spadek temperatury mieszaniny 
gazowej [7, 8] o wartość 4AT5,. 

Po zastąpieniu adiatermicznej osłony diatermiczną 1 po upływie dostate- 
cznie długiego czasu temperatura mieszaniny na skutek wymiany ciepła z 
otoczeniem powróci do wartości początkowej i uzyskany zostanie nowy stan 
równowagi termodynamicznej (ryc. lc). W wyniku izobaryczno-izotermiczne- 
go zmieszania gazów rzeczywistych, objętość powstałej mieszaniny zmieniła 
się o AV w stosunku do sumarycznej objętości jej składników, zaś jej 
entalpia o wartość AH, równą ilości ciepła wymienionego z otoczeniem. Dla 
mieszaniny rzeczywistej równania (2), (4), (6) przyjmą zatem postać 


v(p, R, X) 2 T)+4V(p, ZŚ x); (7) 

h(p, T, x,) = 2% h;(p, T)+4h(p, T, x); (8) 

CAD: JE; X;) a" tp. F)-FAc-(p, A, X;). (9) 
J 


W celu obliczenia poprawki objętości właściwej mieszaniny Av, będącej miarą 
niespełniania się prawa Amagata, posłużono się równaniem stanu gazu 


12 
rzeczywistego w postaci 
RT 
oai; I MOBARDDA (10) 


w którym: B =B(T), C=C(T), D=D(T),... są współczynnikami zależny- 
mi od rodzaju gazu i temperatury. Założono przy tym, że równanie (10) 
ważne jest zarówno dla czystych składników, jak 1 dla mieszaniny. Dla 
mieszaniny współczynniki zależą dodatkowo od jej składu. 


B,„ = B(T, Xx;); C„=C(T, Xx,); D,„ = D(T, Xi); ... (11) 


Po wykorzystaniu zależności (7), (10) i (11) poprawka właściwej objętości 
molowej podczas izobaryczno-izotermicznego mieszania będzie równa 


46 = (B„—) x,B,)+(C„—), x; C;)p+(D„-2 x,D,)p*+ ... (12) 


Znając 40 można w prosty sposób obliczyć poprawkę molowej entalpii 
właściwej Ah. W tym celu posłużono się jedną z podstawowych zależności 


termodynamicznych 
Ch aula z-2P Ę i (13) 
Op T OT p | 


Wykorzystując równanie (7) można wyrażenie (13) doprowadzić do postaci 


0h Óv; (r 
kiej POEZJA RAKA, p pać aso P— 
(6). -Z|s-7(Gr),|-0-7(ar), 


a, kd 
=) x;,|— | +4v-T|-—)|. (14) 
2 (z T OT p 


Całkując równanie (14) wzdłuż izotermy T=idem od ciśnienia p, >0 
do końcowego ciśnienia mieszaniny p, otrzymano 


p OAv 
j 0 0T J, 
Ponieważ przy ciśnieniu p, > O gazy rzeczywiste można traktować jako gazy 
półdoskonałe, zatem 


ho =P x hy. (16) 
J 


p OAv 
h, = 2 X 1+(|40- r(5r),|, dp (17) 
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p Odv 


Wstawiając (12) do (18), otrzyma się po przekształceniach i scałkowaniu 


oraz 


ż | 1 
Ah = [B„— TB„—), x;(B,— TB))] p+z[Cn- TC„-), x;(C;— TC] p*+ 
J J 


l 
+3 LD» TD„-Q x;(D,— TD)]p'+..., (19) 
j 


gdzie 


OB dB, 
= da M B. =-— itd. 
Br e ); dT 


Wykorzystując kolejne równania termodynamiki, wyprowadzono wyraże- 
nie na poprawkę ciepła właściwego, przy stałym ciśnieniu w danej tem- 
peraturze T 


j£ | a. SG > An 
oO ala 06 Wbarz ac 
= — T[(Bn—2 x,B;)p+ż(Cn—2 x,C;)p*+ ...|, (20) 
J J 


gdzie 


| ," 0? Bi, h " d” B, 4 
B,, = b 0T2 0 B; = "m itd. 
jd) 


Postępując podobnie, uzyskano dodatkowo wzory na zmianę objętości 45, 
podczas adiatermicznego zmieszania gazów oraz towarzyszący mu spadek 


temperatury: AT, = T—T, 


HE. a R 
Ava = v„(P. saa X;)-), x; v,(P. T) | e” T)+ 
J 


+8B.(1)-V GB(TH|C.(E)-FoeCAD i "25 
J J 
Ah 
AT,,=T-T, ==G-, (22) 
CHTap 


gdzie Cr, — oznacza średnie molowe ciepło właściwe, przy stałym ciśnie- 
niu w zakresie temperatur od T,, do T. Sposób obliczania średniego ciepła 
właściwego gazu rzeczywistego przy dowolnym ciśnieniu p przedstawiony jest 
w pracach [6, 7]. 
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MIESZANIE GAZÓW RZECZYWISTYCH W WARUNKACH 
IZOCHORYCZNYCH (v = idem) 


Kolejne fazy tworzenia roztworu gazów rzeczywistych przy zachowaniu 
stałej objętości przedstawiono na rycinie 2. W szczelnym i nieodkształcalnym 
zbiorniku, w stanie równowagi termodynamicznej z otoczeniem znajdują się 
składniki gazowe o tych samych ciśnieniach i temperaturach, poprzedzielane 
nieprzepuszczalnymi, diatermicznymi przegrodami (ryc. 2a). Po zaizolowaniu 
cieplnym zbiornika i usunięciu przegród gazy mieszają się ze sobą, tworząc 
roztwór gazowy o parametrach jak na rycinie 2b. Przy dostatecznie dobrej 


a) b) c) 


Ryc. 2 


izolacji termicznej zbiornika i krótkim czasie mieszania się gazów można 
przyjąć, że proces ten przebiega adiatermicznie. Usunięcie izolacji termicznej 
umożliwia wymianę ciepła z otoczeniem i przejście układu do nowego stanu 
równowagi termodynamicznej, co jest jednoznaczne z osiągnięciem początko- 
wej temperatury (przed zmieszaniem), przy jednoczesnym wzroście ciśnienia 
mieszaniny o wartość Ap, właściwej energii wewnętrznej o 4u (ryc. 20). 
Zmiany te spowodowane są również oddziaływaniami międzydrobinowymi. 
W rezultacie, dla mieszaniny gazów rzeczywistych równania (1), (3) i (5) 
przyjmą postać 


P„(V, E, Xj) = p(V,, T)+4p(V, i X;); (23) 

GW, T, x) =Y x;ż;W,, T)+4AG(V, T, x); (24) 
J 

G(W, T, x) = 2, x £,;(V;, T)+A£,(V, T, x,). (25) 
J 


Przyrost ciśnienia Ap = (p„— p)r. ZWany ciśnieniowym efektem izotermiczne- 
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go mieszania, jest miarą niespełniania się prawa Daltona dla gazów rzeczywi- 
stycha ik 2>.żZ]: 

W celu obliczenia Ap, termiczne równanie stanu gazu rzeczywistego (10) 
przekształcono do postaci 


RT Q 
ERN 6 (26) 
D D 
w którym współczynniki wirialne 6, y,ó,... są dla danego gazu zależne 


tylko od temperatury. Dla mieszaniny zależą one dodatkowo od udziałów 
poszczególnych składników. 

Po wykorzystaniu równania (26) ciśnieniowy efekt izotermicznego mieszania 
gazów można obliczyć ze wzoru [7]. 


2 3 
Ap =(Pm-P)r LE 057: „-z|a 5) +Y (5) z e] (27) 


Przy wyprowadzaniu wzoru (27) uwzględniono warunek 


D=0,=),x,0, = idem (28) 


ważny dla szczelnego zbiornika o stałej objętości. 
Zmianę molowej energii wewnętrznej można obliczyć z zależności 


( ALE r( óP 29 
Óujr OT |, zj e 
Podstawiając do (29) ciśnienie i jego pochodną po temperaturze zgodnie z 
równaniem (26) uzyskuje się po odpowiednich przekształceniach i scałkowa- 
niu w granicach od 6, o (p— 0) do v(p, T), następujący wzór na energię 
wewnętrzną gazu 
b PoslPĄ. ME mód JO 
u(p, T) pi 47 —2 + m3 
v(p, T) 25*(p, T) 3v'(p, T) 


w którym apostrofami oznaczono pochodne współczynników wirialnych po 
temperaturze. 
Wstawiając tak wyrażoną energię wewnętrzną do wzoru (24) otrzymuje się 


du(V, p; Xxj) = Um (Dm» I, EDIESLNIA T) <=; 


sM (30) 


Bm= TBm . Ym= Pim 


5 2) Ez Haza Ę 


Zmianę ciepła właściwego przy stałej objętości oblicza się ze wzoru 


4 cAu M p. „6 
4€, = ss; +2 (5 >+ ...||. 32 
Ć - (r), = pz (5 jyga © )- > tz | (32) 


J 


o =P 


z: sez | (31) 


vj J 


Podwójnym apostrofem oznaczono drugie pochodne współczynników wirial|- 
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nych po temperaturze. Temperaturę gazu po adiatermicznym zmieszaniu (ryc. 
2b) wyraża wzór 


(33) 


Cr, 


Wstawiając obliczoną w ten sposób temperaturę mieszaniny do równania 
(26) otrzymuje się ciśnienie w zbiorniku po zmieszaniu gazów rzeczywistych 
bez wymiany ciepła z otoczeniem. 


RT... BAC). MR, 
Da = Ta, Bl Ej SAM 


5 cH Ę% 


WYNIKI OBLICZEŃ I WNIOSKI 


Z wielu postaci równania stanu gazu rzeczywistego wybrano do wstęp- 
nych obliczeń równanie Beattie-Bridgemana 


2 b Cc a 
72 = 0) —- 1-— |-4,(1-= 3 
pv RT|i+B, ( | zs) | s) (35) 


w którym symbolami 4,, B,, a, b, c oznaczono stałe indywidualne rozpatry- 
wanej substancji. 

W literaturze rekomenduje się równanie Beattie-Bridgemana jako jedno z 
lepiej sprawdzonych i dokładnych równań stanu. O wyborze tego równania 
zadecydowała również jego prosta postać oraz przydatność do określania 
stanu mieszanin gazowych. Beattie wykazał doświadczalnie, że zastosowanie 
równania (35) do obliczania objętości właściwej gazów rzeczywistych powo- 
dowało w szerokim zakresie parametrów błąd mniejszy od 0,55'/, [4]. Daje 
się ono również sprowadzić do postaci (10) i (26), wykorzystywanych w 
pracy. W pierwszym przypadku otrzymamy 


pó = RT+B(T)p+C(T)p*+D(T)p', (36) 
gdzie 
B. PJ Pi, KAMMEL ANI 
RER CO= RTZ (KT O) = RT (RT)* " (RT)S 


(37) 


Zależne od temperatury i rodzaju gazu współczynniki wirialne f, y, 0 oblicza 
Się ze wzorów 


-  cR 
5 = —4,+B, RT — za; 
5 B R 
) = aĄ, —bB,RT-—; (38) 
bcB, R 
Ó = o 
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W przypadku mieszanin Beattie zaleca obliczanie stałych 4,, B,, a, b, c ze 
wzorów: 


Apm = [2 Xj (AA: 


kła = > Xj B;;; 
J 


J 
J 

GF =DXCy 
JI 


W drugim przypadku równanie Beattie-Bridgemana przyjmuje postać 


Zaz a OO 
po KTH" (40) 
v U v 


w której współczynniki 5, y, 6 określone są równaniami (38). Wyprowadzone 
równania rozwiązano numerycznie dla różnych składów mieszanin oddecho- 
wych oraz dla parametrów termicznych w zakresie ważności równania 
Beattie-Bridgemana [7, 9]. Niektóre wyniki obliczeń dla najczęściej stoso- 
wanej mieszaniny helowo-tlenowej pokazano na rycinach 3 do 6. Wynika 
z nich, że w zakresie parametrów stosowanych w oceanotechnice, termiczne 
(4p, Ap, AT, AV, AV,, AT.) 1 kaloryczne (Ah, 4a, AC,, 4C,) efekty mie- 
szania przyjmują znaczące wielkości i ich pominięcie na etapie obliczeń 
projektowych może prowadzić do poważnych błędów i przykrych następstw. 

Uzyskane wyniki obliczeń wymagają weryfikacji doświadczalnej. Pod 
względem jakościowym są one jednak zgodne z wynikami badań ekspery- 
mentalnych dla podobnych mieszanin (He — CO»; CO; — NN; N, — H;) 
oraz zbliżone pod względem ilościowym [1, 2]. Niezależnie od badań ekspe- 


2 — Szczecińskie Rocz. Nauk. t. II, z. 1 
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1-P=0750 2.-P-107; 3-P-15-07; 4.-P=2.107 [Pal 
Ryc. 4 


0 02 0,4 06 98 Xe 10 
Ryc. 5 


rymentalnych niezbędne jest również powtórzenie obliczeń za pomocą in- 
nych, dokładniejszych równań stanu. Wyniki takich badań umożliwiają 
opracowanie tablic i wykresów, które ułatwią i uściślą obliczenia cieplne 
hiperbarycznych systemów do badań podwodnych. 
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VUn=V =ldem 
T=273,15K 


Ryc. .6 


WYKAZ OZNACZEŃ 


a 


o B,, «, b, c — stałe indywidualne w równaniu Beattie-Bridgemana 

— molowe ciepło właściwe 

— entalpia 

— molowa entalpia właściwa 

— ciśnienie 

— ciepło 

uniwersalna stała gazowa 

— temperatura absolutna 

— energia wewnętrzna 

— molowa energia wewnętrzna właściwa 

— objętość 

— molowa objętość właściwa 

— udział molowy składnika. 

Indeksy dotyczą: a — procesu adiatermicznego, j — składnika mieszani- 
ny, m — mieszaniny, o — ciśnienia odniesienia, p — stałego ciśnienia, IT — 
stałej temperatury, V — stałej objętości. 
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ROMAN SOBAŃSKI 
THERMAL AND CALORIC EFFECTS OF MIXING REAL GASES 


Summary 


Some problems connected with the mixing of real gases are discussed. It is proved that 
many gas mixtures at parameters observed in practice (pressure < 20 MPa, average temperatu- 
res) should be treated as imperfect solutions. The process of mixing is then accompanied by 


20 


distinct heat effects and variation in thermal parameters. These phenomena occur for instance 
while preparing breathing mixtures used for saturated diving. 

Worked out formulas enable calculation of thermal and caloric effects which accompany the 
process of mixing real gases. Results of calculations for helium-oxygen mixtures quoted in the 
paper, indicate the necessity for taking under account such effects when designing systems for 
underwater exploration. 


POMAH COBAHBCKIH 


TEPMMUECKNE MH KAJIOPHAMETPMUECKNHE J0OQEKTBI 
CMEIIMBAHHMA PEAJIBHbIX I A30B 


Pe3roMe 


B paóoTe OToBopeHbi HeKOTODpbIe IIpOÓJIEMbI, CBA3AHHbIE CO CMELIIHBAHHEM  pEAJIbHbIX 
ra30B. Iloka3aH0, uTo npu rmapaMeTpax, BCTpeuaeMbIX B IIpakTHKE (MaBJIeHHe < 20 Milla, 
yMepeHHbie TeMIIepaTypbl), CMECH MHOTHX TA3OB CJIE[yeT pACCMATpHBATb KAK HEH /1EAJIbHbie 
pacTBOpbi. B 3TOM CJIyuae, Mpouneccy CMEIIIHBAHHA COIIYTCTBYFOT TEIIJIOBbie 3Q)PEKTbI, A TAKE 
H3MeHCHHA TEpMHUECKHX IIADpaMeTpoB. DTM ABJIEHHA HAOJIIONAFOTCA BO BPEMA IIDHTOTOBJIEHHA 
NbIXAaTEJIBHbIX CMECEŃ, ITDHMEHAEMBIX IIDH TJIYOOKOBOJ1HOM HbipAHHH. 

BbIBe1eHbl pacućTHble QOPMYJIBI, IIO3BOJIAKOIIIME IIDOH3BECTH pACUGT TEepMHVYECKHX H KAJIO- 
puMeTpH4ecKkHX 3ć(beKTOB, COIIYTCTBYFOLIHX IIPOLIECCY CMEIIMBAHHA PEAJILHbIX TA30B. IIpuBe- 
AGHHbie B pa6OTe pe3YJIbTATbI PACHETOB JUIA TEJIMO-KACJIODONHbIX CMECeń yKAa3bIBAKOT Ha 
HeO6XOJMMOCTb yućTa 3THX %pbekTOB NpH MpoeKTUDOBAHHH CHCTEM JUIA TJIyYOOKOBOJIHbIX 
HCCJIEHOBAHHA. 
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Streszczenie. W publikacji przedstawiono sposób opracowania wyników 
badań równowag fazowych w wieloskładnikowych układach typu sole — woda. 
Omówiono procedurę wyznaczania najlepszych regresyjnych równań rozpuszczal- 
ności. Wyznaczono równania korelacyjne rozpuszczalności i gęstości roztworów 
nasyconych dla izotermy 293 K układu: FeS0, — NaCl — FeCl, — Na,SO4 — 
HO, na podstawie wyników badań równowagi fazowej w tym układzie. Przedsta- 
wiono sposób wykorzystania równań korelacyjnych rozpuszczalności w rozwiązy- 
waniu zagadnień technologicznych. 


Słowa kluczowe: równowagi fazowe — rozpuszczalność — równania 
korelacyjne układu: FSO, — NaCl — FeCl, — Na,SO, — HO. 


Przemiany fazowe w skondensowanych wieloskładnikowych układach 
typu sole — woda stanowią podstawę wielu operacji technologicznych, 
takich jak rozpuszczanie, krystalizacja, zatężenia roztworów, reakcje podwój- 
nej wymiany itp. Operacje te stosowane są w licznych technologiach prze- 
mysłu nieorganicznego, takich jak przerób surowców mineralnych, otrzymy- 
wanie nawozów mineralnych, otrzymywanie soli i innych związków nieorga- 
nicznych. Prawidłowe zaprojektowanie i prowadzenie procesu w skali prze- 
mysłowej wymaga znajomości równowag fazowych w odpowiednich ukła- 
dach wieloskładnikowych. I tak np. wytwarzanie sody metodą amoniakalną 
oparte jest na wykorzystaniu właściwości fazowych układu: NaCl — 
NH;HCO; — NaHCO; — NH;4CI — H,O. Wytwarzanie nawozów mine- 
ralnych typu NP (fosforany amonowe) wymaga znajomości równowagi fazo- 
wej w układzie: NH; —. P,O; — HO. W przetwórstwie soli mineralnych 
magnezowo-potasowych wykorzystuje się znajomość równowagi fazowej w 
układach: KCI — MgSO4 — K,SO4 — MgCl, — H,O; KCI — MgCl, — 
H,O; K,SO4 — MgSO4 — HO itp. 

Równowagi fazowe w tego typu układach mają często bardzo złożony 
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FeCl-4,0 
N 
08 


06 a 


Ryc. 1. Diagram fazowy układu FeS0, — NaCl — FeCl, — Na,SO4 — H,O w temperaturze 
293 K 


charakter. Jedną z faz jest roztwór nasycony, z którym w równowadze mogą 
występować sole proste i ich wodziany, sole podwójne, roztwory stałe itp. 
Liczba i rodzaj współistniejących w stanie równowagi faz zależy od typu 
układu, składu 1 temperatury. 

Na rycinie 1 przedstawiono wyznaczony doświadczalnie diagram fazowy 
układu: FeS04 — NaCl — FeCl, — Na,SO, — H,O, w temperaturze 
293 K, w układzie współrzędnych Janeckego [1]. W równowadze z roz- 
tworem nasyconym mogą występować następujące fazy: FeSO04:7H,O; 
Na,SO4:10H,O0; Na;SO4; FeSO4:Na,SO4:4H,O; NaCl i FeCl, :4H,O. 
Diagram fazowy stanowi jedynie graficzną ilustrację złożonego charakteru 
fazowego układu. Przedstawiono na nim przebieg linii rozgraniczających 
pola krystalizacji poszczególnych faz (linie eutoniczne) i punkty zerozmienne. 
W omawianym przypadku stanom dwufazowym układu (roztwór nasycony 
+ faza stała) odpowiadają w trójwymiarowej przestrzeni składu obszary w 
kształcie ostrosłupów, ograniczonych powierzchniami składu roztworów na- 
syconych (tzw. pola krystalizacji), których wierzchołki leżą w punktach 
składu krystalizujących soli. W ustalonych warunkach ciśnienia i temperatu- 
ry maksymalna liczba mogących współistnieć w stanie równowagi faz wynosi 
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cztery, co odpowiada punktom zerozmiennym E,—Ey, zaznaczonym na 
diagramie fazowym (ryc. 1). W miarę wzrostu liczby składników w układzie, 
jego właściwości fazowe stają się coraz bardziej złożone, a ich odwzorowanie 
w formie diagramów fazowych jest coraz trudniejsze. 

Projektowanie procesu technologicznego wymaga znajomości stanu fazo- 
wego układu dla różnych stężeń składników i temperatury. Stan fazowy 
układu określony jest wtedy, gdy dla danych warunków stężeń składników i 
temperatury można określić liczbę i skład współistniejących w stanie równo- 
wagi faz, a zwłaszcza skład roztworu nasyconego. 

Obliczenia technologiczne procesów odparowania, krystalizacji, wymiany 
jonowej itp., wykonuje się najczęściej w sposób przybliżony, metodami 
graficznymi, wykorzystując wykresy pomocnicze zależności składu roztworu 
nasyconego od stężeń poszczególnych składników. W publikacji przedstawio- 
no sposób opracowania wyników badań równowag fazowych w układach 
wieloskładnikowych i wyznaczania równań korelacyjnych rozpuszczalności 
metodami analizy regresji. Tego typu równania pozwalają na znaczne 
usprawnienie obliczeń technologicznych procesów związanych z przemia- 
nami fazowymi. Proponowaną metodę zilustrowano na przykładzie izotermy 
rozpuszczalności czteroskładnikowego układu: FeS04 — NaCl — FeCl, — 
Na,;SO4 — H,O. Ponieważ metoda wyznaczania równań korelacyjnych 
rozpuszczalności może być zastosowana do układów o dowolnej liczbie 
składników i w dowolnych układach współrzędnych stosowanych przy spo- 
rządzaniu diagramów fazowych, dla jednoznacznego przedstawienia proble- 
mu konieczne jest wprowadzenie pewnego formalizmu matematycznego w 
zakresie algebry przestrzeni liniowych 1 geometrii analitycznej przestrzeni 
wielowymiarowych [2]. 


Sformułowanie problemu. Skład układu s-składnikowego może być 
odwzorowany w (s-l) wymiarowej przestrzeni x, zwanej dalej przestrzenią 
składu. Współrzędne dowolnego punktu 


x(X4, X2 wee j X 11) 


tej przestrzeni określają stężenia poszczególnych składników w układzie. 
Przestrzeń liniowa X rozdzielona jest na dwie części, odpowiadające stanom 
jednofazowym i wielofazowym układu. Właściwości fazowe układu są okreś- 
lone wtedy, gdy każdemu punktowi x,, leżącemu w obszarze wielofazowym 
przestrzeni składu, można przyporządkować jednoznacznie punkty X, i Xp, 
określające składy pozostających w równowadze faz. W układach skondenso- 
wanych jedną z faz jest roztwór nasycony. Punkty składu roztworów nasyco- 
nych leżą na pewnej (s—2) — wymiarowej hiperpowierzchni, w (s—1) — 
wymiarowej przestrzeni składu. Hiperpowierzchnia ta jest miejscami gładka 
w obszarach odpowiadających krystalizacji poszczególnych faz stałych $,. Na 
hiperpowierzchni składu roztworów nasyconych stężenia składników powią- 
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zane są pewną zależnością funkcyjną F (xp) = O, O nie znanej na ogół postaci 
analitycznej. 

Na rycinie 2 przedstawiono tego typu odwzorowanie właściwości fazo- 
wych układu czteroskładnikowego w trójwymiarowej przestrzeni składu. 
Może to być np. układ typu: AX — BX — CX — H;,O, gdzie symbolami 


Ryc. 2. Diagram fazowy czteroskładnikowego układu typu AX — BX — CX — HO 


AX, BX, CX — oznaczono sole o wspólnym anionie. Na osiach układu 
współrzędnych odłożone są stężenia trzech dowolnych składników układu. 
Powierzchnie H;,, H,, Ha określają skład roztworów nasyconych względem 
faz S,, S, i Sz. Punkty E;>, E3, E,3 określają skład roztworów nasyconych 
dwiema solami w granicznych układach trójskładnikowych (punkty eutoni- 
czne). Punkt E,;s jest punktem zerozmiennym w danej temperaturze i 
określa skład roztworu nasyconego trzema solami. Odcinek R(x,, x,) wyzna- 
czony przez punkty x, i x, nosi nazwę promienia krystalizacji. Jeżeli będą 
znane równania powierzchni składu roztworów nasyconych H;, H, i Ha, dla 
dowolnego składu układu w obszarze krystalizacji fazy S,, skład roztworu 
nasyconego xp moŻNa wyznaczyć znajdując punkt przebicia powierzchni H, 
przez promień krystalizacji R(x,, x,). Układ współrzędnych przestrzeni skła- 
du zależy od sposobu wyrażania stężeń. W przypadku układów typu sole — 
woda, w których sole mają różne aniony, powszechnie przyjętym sposobem 
wyrażania stężeń jest tzw. układ współrzędnych Jiineckego, w którym skład 
układu wyrażany jest przez ułamki molowe jonów występujących w układzie 
i liczbę molową wody. Ten układ współrzędnych jest szczególnie wygodny 
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dla wyznaczania równań korelacyjnych rozpuszczalności, gdyż powierzchnie 
składu roztworu nasyconego można aproksymować równaniami postaci: 


M = f(x, Xą, ::., Xs— 2) (1) 


gdzie wektor (x,, Xa, ..., X;_2) Określa ułamki molowe jonów w roztworze 
nasyconym, natomiast M oznacza liczbę molową wody. 

W analogiczny sposób można przy tym skorelować właściwości fizyko- 
chemiczne roztworu jak np. gęstość, ze stężeniami soli w roztworze nasyco- 
nym: 


d=g(x,, Xą, ..:, Xs- 2): (2) 


W ten sposób równanie 1 określa skład roztworu nasyconego jako funkcję 
stężeń rozpuszczonych soli, natomiast równanie 2 określa zależność gęstości 
roztworów nasyconych od tych samych zmiennych niezależnych [3]. 


Wyznaczanie równań korelacyjnych. Rozpuszczalność soli i rów- 
nowagi fazowe w układach wieloskładnikowych wyznacza się metodami 
doświadczalnymi. W tym celu dla różnych punktów składu układu wyznacza 
się skład roztworu nasyconego, często również jego gęstość i skład pozostają- 
cej z roztworem w równowadze fazy stałej. W ten sposób wyznaczone zostają 
wartości zmiennych (M,, d,) dla różnych punktów składu określonych wekto- 
rami (X4, Xą, ..., Xs_2) 1 dla obszaru krystalizacji fazy $,. Ponieważ postać 
analityczna zależności funkcyjnych określonych równaniami 1 i 2 nie jest 
znana, można ją aproksymować dowolnymi funkcjami analitycznymi. Naj- 
prostsze rozwiązanie polega na aproksymowaniu tych równań wielomianem 
(s—2) zmiennych niezależnych: 


M, =),ayXq X7...X$_2, (3) 
l 
— r S t 
PPI EAĄCZZNECA CE (4) 
k 
gdzie: m, n, p, r, s, t =O, 1, 2,...; natomiast indeks j — oznacza, że równanie 


aproksymujące odnosi się do obszaru krystalizacji fazy S,. Występujące w 
tych szeregach współczynniki regresji a, 1 by można wyznaczyć znając 
wartości M, i d,,, gdzie i — oznacza numer pomiaru. Jeżeli przyjmie się, że 
zmienne M i d są niezależne, wszystkie M, i d,, mają jednakowe wariancje i 
pochodzą z populacji o rozkładach normalnych, natomiast macierz (X: X7) 
jest nieosobliwa, to nieobciążonymi i najefektywniejszymi estymatorami 
współczynników regresji będą wektory 4 i b, wyznaczone metodą najmniej- 
szych kwadratów: 


i=(X:X7) ':X:M (5) 
b=(X:X7) "-X:d (6) 
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gdzie: M =(M;, M,,..., M,)'; d=(dy,d,,...,d,)”;n — liczba pomiarów; 
X — macierz pomiarów o wymiarach n x(s— 2). 
Oczywiste jest, że równania regresji należy wyznaczyć dla obszarów krystali- 
zacji wszystkich faz. Pozostaje do rozstrzygnięcia stopień wielomianu w 
równaniu regresji i liczba uwzględnionych w równaniu wyrazów. 
Optymalną postać równania regresji można określić testując istotność 
poszczególnych wyrazów w równaniach (5, 6). W tym celu jako pierwsze 
przybliżenie przyjmuje się pełne wielomiany pierwszego, drugiego, trzeciego 
itd. stopnia (s— 2) zmiennych niezależnych. Istotność każdego z wielomianów 
testuje się na podstawie stosunków średnich kwadratów odchyleń zniesio- 
nych przez wyrazy danego stopnia do wariancji resztowej wyznaczonego 
równania i wariancji metody doświadczalnej (wariancyjny test istotności 
korelacji). Po określeniu optymalnego stopnia wielomianu testuje się istot- 
ność każdego wyrazu w równaniu w oparciu o statystykę 


2 
(= L 
Sbl 
gdzie: b, — l-ty współczynnik regresji, sbi wariancja współczynnika 
regresji. Statystyka t* ma rozkład Fishera — Snedecora o 1 i n-l—1 
stopniach swobody (I — liczba estymowanych współczynników regresji). 


Jeżeli nie jest spełniona relacja t;? > F, y „_,-4, wÓwczas odrzuca się testowa- 
ny wyraz z równania regresji. Taka procedura pozwala na ograniczenie 
efektów nadmiernego dopasowania równań regresji, co umożliwia ich wyko- 
rzystanie do ekstrapolacji danych doświadczalnych. 


Omówienie wyników. Przedstawioną procedurę wyznaczania równań 
rozpuszczalności soli 1 gęstości roztworów nasyconych zastosowano do opra- 
cowania wyników badań równowagi fazowej w układzie FeS0, — NaCl — 
FeC!1, — Na,S04 — H,O w temperaturze 293 K. W tabeli 1 zestawiono 
wyniki sekwencyjnej analizy regresji dla rozpuszczalności soli w obszarze 
krystalizacji fazy FeSO,4:Na,SO04:4H,O, a także wartości współczynnika 
korelacji wielokrotnej R, wariancji resztowej s” i wartości testu Fishera— 
Snedecora dla testowanego równania. Z wartości tam zestawionych wynika, 
że istotny wkład w równania korelacyjne wnoszą wyrazy pierwszego 1 
drugiego stopnia (współczynniki regresji ag — a, Oraz aa — as). Spośród 
wyrazów trzeciego stopnia najistotniejsze są wyrazy a;y” i agx*, natomiast 
pozostałe wyrazy należy odrzucić jako nieistotne. Analizą objęto 11 równań 
regresji, spośród których za najlepsze należy uznać równanie oznaczone 
numerem 9. Charakteryzuje się ono wysoką istotnością, na co wskazuje 
wartość współczynnika korelacji wielokrotnej (R* = 0,90901) i niewielka wa- 
riancja resztowa. Wartość testu F dla tego równania wynosi 288,05, przy 
wartości granicznej F dla x = (0,05 wynoszącej 2,16. 


2% 


Tabela 1. Analiza równań regresji dla obszaru krystalizacji soli 
FeSO4 : Na,SO, *4H,O w temperaturze 293 K 


Ogólna postać wielomianu: 
M =aqg+a, x+tazy+azx*+a4y?+asxy+agx +ay*+agx* y+dag xy? 


gdzie: M — liczba molowa wody; x — ułamek molowy Fe**; y — ułamek 
molowy CI 


Testowane Liczba 
współczynniki stopni 
regresji swobody 


A+, d>/aQ 0,37469 0,59617 
az — as/ay — A 0,88380 ATZ 
A> — dg/dg — A> 0,89088 0,10646 
A> —do/dg —aG 0,90905 0,09032 


dułks=d< 0,88949 0,10782 
ZA 0,89823 0,09929 
g/g A> 0,88431 0,11288 
49/ag — 3 0,88529 0.1! 130 


Aa dx, 2 0,90901 0,08930 
Ge faq dt.A> 0,89944 0,09869 
ME YA 0,89823 0,09988 


l 
2 
3 
4 
5 
6 
7 
8 
9 
0 
1 


— M 


Należy podkreślić, że przedstawiona procedura wyznaczania najlepszego 
równania regresji pozwala na przetestowanie szeregu równań, bez konie- 
czności wielokrotnego powtarzania obliczeń. W tym celu macierz (X:X7) 
oblicza się dla dostatecznie wysokiego stopnia wielomianu (zazwyczaj trzecie- 
go lub czwartego), a następnie eliminuje się wyrazy statystycznie nieistotne. 
W podobny sposób wyznaczono równania rozpuszczalności dla obszarów 
krystalizacji wszystkich faz, które ilustruje tabela 2. Z, kolei w tabeli 3 
zestawiono równania zależności gęstości roztworów nasyconych od ich skła- 
du, będące modelem regresyjnym izotermy rozpuszczalności soli w układzie: 
FeS04 — NaCI — FeCl, — Na,SO4 — H,O, w temperaturze 293 K. 

Wyznaczone równania można wykorzystać w modelowaniu procesów 
odparowania roztworów, krystalizacji soli, konwersji siarczanu żelazawego z 
chlorkiem sodowym itp. Na rycinach 2 i 3 przedstawiono rzuty poziome 
izotermy rozpuszczalności układu: FeS04 — NaCl — FeCl, — Na;SO4 — 
H,O, z naniesionymi liniami stałych wartości liczby molowej wody (izohyd- 
ry) i gęstości roztworów nasyconych (izopykny), wykreślonymi z wykorzysta- 
niem równań regresji. Linie rozgraniczające poszczególne obszary krystaliza- 
cji zaznaczone na tych wykresach zostały wyznaczone jako linie przecięcia 
powierzchni składu roztworów nasyconych. Są one rozwiązaniami układów 
równań zestawionych w tabeli 2. 

Układ czteroskładnikowy: FeS0, — NaCl — FeCl, — Na,;SO4 — H,O, 
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Tabela 2. Zestawienie współczynników regresji równań korelacyjnych rozpuszczalności 
dla izotermy 293 K układu FeS0, — NaCl — FeCl, — Na,SO, — H,O 


M =ag+a;,x+a,y+az x +aąy? tas xy+agx +ay +agx*y+dgxy*, 


gdzie: M — liczba molowa wody; x — ułamek molowy Fe**; y — ułamek molowy CI 


l 
2 
3 
4 
5 
6 
7 
8 
9 
0 


l 


Współ- 
czynnik 
regresji 


14,58565 
19,26639 
5,12303 

— 30,98964 
5,39616 

— 6,47811 
16,95478 
— 10,57356 


Obszary krystalizacji soli 


ER AGU eSO,:7H,O 9050,10; NI |Feiylo 
PAA 


—0,17808 
19,22974 

— 3,07211 
— 3,4801 1 
27,22924 
— 1,95817 
—32,12952 


— 3,00925 


20,70219 
=23,261.30 
— 22,56161 

1,03614 

65,91249 

26,05662 


— 80,08858 


—65,09607 


5,38405 
—0,91707 
3,67715 
1,06906 
=©52271 


—4,08577 


—(0,15895 


— 213453 
—18,73136 | 
27,61087 
—0,30032 
— 24,25274 
23,52331 


— 
— 


liczba 
pomiarów 120 94 74 


0,90901 0,90041 0,96857 0,87276 
0,08930 1,06149 0,36867 0,17841 
288,05 144,65 378,60 76,59 


0,96769 
0,00178 
101,83 


jest tak zwanym układem wymiennym. Możliwe jest w nim zachodzenie 
reakcji podwójnej wymiany zgodnie z równaniem: 


FeSO, + 2N aCl =FeCl, + Na,SO4. 


Z diagramu fazowego układu wynika, że w temperaturze 293 K tworzący się 
w wyniku konwersji siarczan sodowy współkrystalizuje z siarczanem żelaza- 
wym jako sól podwójna: FeS$O4 : Na,SO, :4H>0O. Z technologicznego punktu 
widzenia interesująca jest wydajność reakcji konwersji względem jonu sodo- 
wego W,,+. Wydajność ta zależy w danej temperaturze od stosunku molo- 
wego FeS$O04 do NaCl oraz ilości wody wprowadzonej do układu. Zależność 
tę można określić wykorzystując wyznaczone równania korelacyjne rozpu- 
szczalności. 

Na rycinie 4 przedstawiono komputerowy wykres zależności Wy,+ od 
ułamka molowego FeSO, i liczby molowej wody w mieszaninie reakcyj- 
nej. Z wykresu wynika, że maksymalną wydajność soli podwójnej 
FeSO4 : Na,SO, *4H,O, wynoszącą 0,65— 0,70, można uzyskać z mieszaniny 
reakcyjnej zawierającej 0,6 mola FeSO4, 0,4 mola NaCl i 3,5—4 moli wody. 
W tych warunkach 65—70"/, chlorku sodowego ulegnie konwersji do soli 
FeSO4 : Na,SO4 *4H,O, która wykrystalizuje z roztworu. 
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Tabela 3. Zestawienie współczynników regresji równań korelacyjnych gęstości roztworów nasy- 
conych dla izotermy 293 K, w układzie FeS0, — NaCl — FeCl, — Na;SO4 — H,O 


d = bo +b, x +b, y+bz x? +b, y? +bs xy +bgx* +b> p? +bg x” y+bg xy? 


gdzie: d — gęstość roztworu nasyconego kg/dm*; x — ułamek molowy Fe**; y — ułamek 
molowy CI" 


Współ- Obszar krystalizacji soli 


czynnik , j 
cares odj FeSO,-7H,O|Na„SO,-10H,O FeCl,-4H,O 
2 


1,340386 1,60001 1,15122 1,85789 2,02691 
—0,182933 —0,61736 0,12636 0,37681 — 1,78246 
—0,026323 0,027138 —0,05158 — 1,26350 0,98224 

0,505202 0,255210 0,52820 —0,25449 0,78249 
— 0;279694 — 1,041420 0,14363 0,60510 —0,83434 
— 0,355942 0,150290 0,10478 —0,20358 0,26986 
—0,133151 = = 0,39660 — 

0,246683 1,19316 0,16921 
—(0,000027 a =] 

0,457693 nĘ = 


© O JOU A LN — 


liczba 
pomiarów 139 115 99 


R* 0,96228 0,87568 0,96553 0,99584 0,98950 
s” 10% 0,58 1,10 0,30 0,03 
Fot 365,66 320,41 


yz 


Fe 90 Na* Na„ SO, 


Ryc. 3. Izoterma 293 K układu FeSO, — NaCl — FeCl, — Na,SO4 — H,O z naniesionymi 
izohydrami 
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FeCl, ge | | R 


0,8 


Kor 
04 4 
Ś 
0,2 
9 
FeS0y Na50, 


Ryc. 4. Izoterma 293 K układu FeSO, — NaCl — FeCl, — Na;SO04 — HO z naniesionymi 
izopyknami 


A 
FESO, 


Ryc. 5. Zależność wydajności Wy,+ reakcji podwójnej wymiany FeSO, z NaCl od składu 
mieszaniny reakcyjnej 
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Przedstawione przykłady nie wyczerpują wszystkich możliwości wykorzy- 
stania wyznaczonych równań korelacyjnych. Zaletą ukazanej tu metody 
opracowania wyników badań równowag fazowych (rozpuszczalności) w ukła- 
dach wieloskładnikowych jest fakt, że na tym samym zbiorze punktów do- 
świadczalnych można korelować według jednolitej procedury zarówno rozpu- 
szczalności soli, jak i dane dotyczące właściwości fizykochemicznych roztwo- 
rów nasyconych. Sekwencyjna analiza istotności wyrazów w równaniu re- 
gresji pozwala przy tym na ograniczenie efektów nadmiernego dopasowania 
równań do zbioru danych. 

Wyznaczone równania stanowią zwięzłe ujęcie ilościowe wyników badań 
fizykochemicznych i pozwalają na znacznie efektywniejsze ich wyzyskanie do 
rozwiązywania zagadnień technologicznych. 


ALEKSANDER PRZEPIERA, STANISŁAW BURSA, 
KRYSTYNA PRZEPIERA 


CORRELATION OF SOLUBILITY IN MULTI-COMPONENTAL 
SYSTEMS OF THE TYPE SALTS — WATER 


Summary 


The way of describing phase equilibrium results in multi —componental systems of the type 
salts-water is presented in the paper. The paper discusses the procedure of calculating the best 
regression equations of solubility. A correlating equation of solubility and density of saturated 
solutions for the isotherm 293 K in the system FeSO0, — NaCl — FeCl, — Na;SO4 — H,O is 
calculated. The calculation is based on the results of research conducted on the phase 
equilibrium in this system. The way of employing correlative equations of solubility in working 
out technological issues is presented. 


AJIEKCAHĄOP IIIOIIEPA, CTAHHCJIAB BYPCA, KPBICTBIHA INIHIDIIEPA 


KOPPEJIAIIMA PACTBOPAEMOCTM B MHOTOKOMIIOHEHTHEBIX 
CHCTEMAX TUIIA COJIM — BOA 


Pe3roMme 


B ny6rnnkanuu IpercTaBJieH Ccroco6 o6paO0oTkH pe3YJIBTATOB MCCJIEHTOBAHHŃ (bA3OBBIX 
paBHOBECcHi B MHOTOKOMIIOHEHTHbIX CHCTEeMaX THNA COJIM — BO.IA. OroBopeHa npoledypa 
BbIBEJTICHHA HAHJIYULNIMX PerpeCcHOHHbIX yIIDABHCHHA PACTBOPAEMOCTH. BbIBeTeHbl YpABHCHHA 
KOPpEJIAIHM DPACTBODAEMOCTH MH NJIOTHOCTH HACbHNEHHbIX PpACTBODOB JIA MZOTEpMbi 293 K 
CHCTEMbl FeSO0, — NaCl — FeCl, — Na,SO, — H,„O, Ha OcHOBe pe3yJIBTATOB MCCJIETOBAHHŃ 
(ba30BOTO paBHOBECHA B 3TOŃ cHcTeMme. IlpencTaBJIeH CIIOCO6 MKCIOJHbZOBAHHA YpABHEHHA KOpP- 
peJIdiiMM pAaCTBOpAEMOCTM IIDH DEIIEHMM TEXHOJIOTHYECKHX IIPDOÓJIEM. 
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FRYDERYK STRĘK 


BADANIA DOŚWIADCZALNE WNIKANIA CIEPŁA 
W MIESZALNIKACH CIECZY 


Instytut Inżynierii Chemicznej i Chemii Fizycznej Politechniki Szczecińskiej, 
Al. Piastów 42, 71-065 Szczecin 


Streszczenie. W artykule przedstawiono rezultaty badań własnych wnika- 
nia ciepła w mieszalnikach cieczy dla różnych typów mieszadeł szybkoobrotowych. 
Pomiary wykonano dla burzliwego zakresu przepływu cieczy w zbiorniku mie- 
szalnika. Do badań użyto różnych cieczy niutonowskich oraz czterech zbiorników 
o średnicach od 0,15 m do 0,9 m, z dnem płaskim oraz przegrodami. 

Opierając się na wynikach równoczesnych pomiarów wymiany ciepła oraz 
mocy mieszania dokonano oceny efektywności mieszania różnych mieszalników 
oraz zaproponowano rozwiązania najlepsze dla takiego przypadku mieszania. 


Słowa kluczowe: mieszanie cieczy — mieszadła — mieszalniki — wnikanie 
ciepła — efektywność mieszania cieczy. 


WSTĘP 


Aparaty typu zbiornikowego zaopatrzone w mieszadło mechaniczne, zwa- 
ne ogólnie mieszalnikami, stosowane są w przemyśle chemicznym i prze- 
mysłach pokrewnych do realizacji różnych procesów technologicznych. Zwy- 
kle w praktyce przemysłowej aparaty takie mają nazwy bardziej szczegółowe, 
np. reaktory, lub też w przypadku stosowania ich do konkretnych procesów: 
sulfonatory, nitratory, krystalizatory, absorbery itp. 

Projektowanie mieszalników napotyka do dnia dzisiejszego trudności 
ze względu na złożoność procesów zachodzących w takich aparatach. Można 
zaryzykować twierdzenie, że dokładność obliczeń projektowych aparatów z 
mieszadłami jest mniejsza niż innych aparatów stosowanych w przemyśle 
chemicznym. Również powiększanie skali takich aparatów nie jest dokładnie 
rozeznane szczególnie dla procesów z udziałem reakcji chemicznej. Przyczyna 
takiego stanu rzeczy tkwi z jednej strony w złożonej hydrodynamice miesza- 
nia cieczy w takich aparatach, z drugiej zaś strony w dużej liczbie odmian 
konstrukcyjnych mieszalników, co utrudnia ich gruntowne przebadanie oraz 
optymalizację. Jako dodatkową trudność przy ocenie pracy mieszalników 
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oraz ich projektowaniu można wymienić małą liczbę kompletnych badań 
takich aparatów prowadzonych pod kątem rozwiązań optymalnych [11], co 
powoduje konieczność powtarzania od nowa pewnych fragmentów wykona- 
nych już przez innych autorów pomiarów. 

Wymiana ciepła w mieszalniku jest ważnym procesem decydującym nie- 
kiedy o prawidłowej pracy aparatu. Dotyczy to tych przypadków, kiedy 
mieszalnik jest reaktorem chemicznym 1 zachodzi konieczność odprowadze- 
nia z układu znacznego ciepła reakcji. Zapewnienie odpowiedniej intensyw- 
ności wymiany ciepła w mieszalniku limituje więc wydajność technologiczną 
aparatu. 

Intensyfikacja odbioru ciepła od cieczy mieszanej przy ustalonych tempe- 
raturach skrajnych, oraz ustalonej powierzchni wymiany ciepła, sprowadza 
się do intensyfikacji wnikania ciepła pomiędzy mieszaną cieczą oraz Ściankę 
aparatu lub wężownicy. Znajomość równań opisujących współczynnik wni- 
kania ciepła w szerokim zakresie zmian parametrów geometrycznych układu 
oraz fizycznych i kinematycznych procesu jest więc bardzo ważna. Opubli- 
kowane na ten temat informacje [3, 4, 6] zawierają luki i nieścisłości, co było 
między innymi przyczyną podjęcia na ten temat własnych badań doświad- 
czalnych. 


CELE I ZAKRES BADAŃ WŁASNYCH 


Badania doświadczalne wnikania ciepła w mieszalnikach cieczy prowa- 
dzone są w Zakładzie Inżynierii i Aparatury Chemicznej Politechniki Szcze- 
cińskiej od wielu lat [2, 5, 7—10, 12, 15, 16]. 

Cele badań są wielorakie, a mianowicie: 

1) zaproponowanie — na podstawie własnego materiału doświadczalnego 
możliwie dokładnych do celów projektowych równań do obliczeń współ- 
czynnika wnikania ciepła w mieszalnikach cieczy w szerokim zakresie zmian 
parametrów procesu; 

2) wyciągnięcie — na podstawie własnych wyników badań — wniosków co 
do optymalnych rozwiązań mieszalników dla wnikania ciepła; 

3) testowanie nowych konstrukcji mieszalników zaproponowanych przez 
innych autorów bądź też propozycji własnych. 

Badania wykonano dla zbiorników z dnem płaskim oraz z przegrodami, 
ogrzewanych na części cylindrycznej zaopatrzonych w mieszadła szybko- 
obrotowe (ryc. 1). Pomiary przeprowadzono dla burzliwego zakresu prze- 
pływu cieczy w zbiorniku mieszalnika. Stosowano różne średnice mieszalni- 
ków: D=0,15 m; 0,3 m; 0,45 m; oraz 0,9 m; różne typy mieszadeł 
szybkoobrotowych: turbinowe tarczowe, turbinowe o pełnej długości łopatek 
o różnym pochyleniu oraz mieszadła śmigłowe (ryc. 2); również różne ciecze 
mieszane: woda wodociągowa, woda destylowana, roztwory wodne gliceryny 
oraz różne oleje maszynowe. 
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Ryc. 2. Podstawowe wymiary przebadanych mieszadeł: a) mieszadło turbinowe tarczowe, 
b) mieszadło turbinowe o łopatkach pochylonych, c) mieszadło śmigłowe 


Pomiary prowadzono zarówno w seriach tradycyjnych jak też dla pew- 
nych mieszadeł jako eksperymenty planowane. 


OPIS APARATURY DOŚWIADCZALNEJ ORAZ SPOSOBU 
PROWADZENIA POMIARÓW 


Podstawową częścią zestawu doświadczalnego (ryc. 3) był zbiornik mie- 
szalnika (1) z płaszczem grzejnym (2) i przegrodami (3), zaopatrzony w 
mieszadło szybkoobrotowe (4). Mieszadło napędzane było silnikiem prądu 
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odpowietrzenie 
|. — e: —e: > 
płaszcza 


kondensat 


Ryc. 3. Schemat aparatury doświadczalnej: ! — zbiornik mieszalnika, 2 — płaszcz grzejny, 
3 — przegrody, 4 — wał napędowy z mieszadłem, 5 — silnik napędowy prądu stałego, 
6 — urządzenie sterownicze i pomiarowe mocy oraz częstości obrotów, 7 — zbieracz prądu, 
8 — tuleja tensometryczna, 9 — pompa obiegowa, 10 — zbiornik pośredni cieczy mieszanej, 


11 — chłodnica, 12 — zespół termoelementów, 13 — rejestratory, 14 — garnki kondensacyjne 


stałego (5). Płynną regulację częstości obrotów silnika zapewniał zespół 
napędowy (6), zbudowany w Instytucie Elektrotechniki Politechniki Szcze- 
cińskiej, liczbę obrotów mieszadła zliczano elektronicznie, stosując jako 
czujniki impulsów kontaktrony. Do pomiaru momentu skręcającego mie- 
szadła wypróbowano dwie metody: elektryczną oraz tensometryczną. Metoda 
tensometryczna okazała się bardziej dokładna i dlatego tę metodę stosowano 
w większości wykonanych pomiarów. Do badań wykorzystano zbudowany 
we własnym zakresie torsjometr tensometryczny (8), przystosowany do po- 
miaru małych momentów skręcających oraz chroniący tuleję tensometryczną 
przed siłami bocznymi działającymi na wał mieszadła. 

Pomiary współczynnika wnikania ciepła prowadzono w warunkach usta- 
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lonej wymiany ciepła. Warunki te zapewniano stosując zewnętrzny obieg 
cieczy mieszanej. Ciecz ta spływała z mieszalnika (1) przelewem do zbiornika 
pośredniego (10), skąd zasysana była pompą (9) oraz przetłaczana poprzez 
chłodnicę (11) i rotametr ponownie do mieszalnika (1). Strumień cieczy 
będącej w obiegu musiał być dostatecznie mały w porównaniu ze strumie- 
niem cyrkulacji wewnętrznej cieczy w mieszalniku, wywoływanym mieszad- 
łem, aby nie spowodować wyrażnego zaburzenia pola temperatur cieczy w 
aparacie, a tym samym wymiany ciepła. Wielkość tego strumienia cieczy 
ustalono eksperymentalnie dla każdej wielkości mieszalnika. Wynosił on 
około '/29 do 1/59 strumienia cyrkulacji wewnętrznej. 

Temperatury cieczy: na wlocie do mieszalnika i na jego wylocie, jak 
również średnią temperaturę cieczy w zbiorniku mieszalnika oraz pary 
grzejnej mierzono za pomocą termoelementów (12) podłączonych do rejestra- 
tora (13) oraz kontrolowano za pomocą termometrów rtęciowych. Do ogrze- 
wania mieszalnika stosowano nasyconą parę wodną, wytwarzaną w kotle 
parowym opalanym olejem napędowym. 


OBLICZANIE WSPÓŁCZYNNIKA WNIKANIA CIEPŁA 
NA PODSTAWIE DANYCH POMIAROWYCH 


Współczynnik wnikania ciepła a obliczano dla każdego pomiaru z równa- 
nia definicyjnego współczynnika przenikania ciepła k 


TIGAR (0 
k 4 A, 
gdzie: 
Q 
k ==; 2 
PG "G,) (2) 


Strumień ciepła dostarczany do cieczy mieszanej wynosił 
Q=Gc(t, —t4). (3) 


Współczynnik wnikania a, dla kondensującej się pary na pionowej po- 
wierzchni płaszcza obliczano z równania 


K 
3 'H R q 
gdzie: stała K = f(t,) jest funkcją temperatury nasycenia pary t,. 
Równanie (4) sprawdzono mierząc temperatury Ścianki wewnętrznej mieszal- 
nika i stwierdzono dobrą jego zgodność z doświadczeniem [5]. Uściślona 


w ten sposób stała K w stosunku do danych literaturowych [2] dla zakresu 
t, =(100—103)?C wynosi: K =66101 ty**. 


(4) 


Oy = 
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W obliczeniach wykorzystano specjalnie uprzednio wykonane charaktery- 
styki natężenia przepływu G cieczy przez rotametry oraz charakterystyki 
temperaturowe termoelementów do wyznaczenia napędowej różnicy tempera- 
tur (1,—t,„) przenikania ciepła. 


UOGÓLNIENIE WYNIKÓW BADAŃ WSPÓŁCZYNNIKA 
WNIKANIA CIEPŁA 


Współczynnik wnikania ciepła w mieszalnikach cieczy opisuje się zwykle 
równaniem kryterialnym 


Nu = CRe* Pr" V* (5) 


słusznym dla rodziny mieszalników geometrycznie podobnych o stałych 
wartościach inwariantów podobieństwa geometrycznego. Odstępstwo od tej 
geometrii parametrów a/d, b/d, h/d, ..., powoduje zmianę stałej Ć, co można 
uwzględnić wprowadzając do tego równania poprawkę Y = Y,Y,Y,„... Jeśli 
przyjmiemy określone proporcje wymiarów liniowych mieszalnika za stan- 
dardowe , = (3); . = (al - = (3): ..., dla których będzie € =C,, otrzy- 
mamy dla dowolnej innej geometrii równanie 


Nu = C,Re* PrEVFY. (6) 


Liczby kryterialne występujące w równaniach (5, 6) mają dla mieszalni- 
ków następujące definicje: 

2 
NpEL RAEC Pr =— Oraz V=—. (7) 
A 7 A Hs 

Liczba Reynoldsa została zdefiniowana na podstawie prędkości obwodo- 
wej mieszadła (ndn) oraz wymiaru liniowego (d), czyli ma znaczenie zastępcze, 
gdyż nie jest związana z rzeczywistym polem prędkości cieczy w mieszalniku. 
Różne mieszadła o tych samych średnicach oraz pracujące z tą samą 
częstością obrotów wytwarzają w tym samym zbiorniku i tej samej cieczy 
różne pola prędkości, co powoduje, że równanie (6) nie jest uniwersalne. 
Stałą C, oraz poprawkę Y należy więc wyznaczać indywidualnie dla każdego 
układu mieszadło — zbiornik. 

Z tej przyczyny wielu autorów proponuje inne definicje liczby Reynoldsa 
dla mieszalników, lepiej opisujące burzliwość cieczy mieszanej. Jedną z cie- 
kawszych propozycji, która znalazła uznanie, jest wprowadzenie do tej 
definicji energii jednostkowej mieszania (P/V), czyli energii rozpraszanej przez 
mieszadło w jednostce objętości cieczy. 

Nową definicję otrzymamy [3, 6] podstawiając za częstość obrotów (n) 
wyrażenie powstałe z przekształcenia liczby mocy Ne =P/n*d*og. Przyjmuje- 
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my ponadto: Ne = const, d/D =const, H/D = const, oraz V =nD*”H/4, 
wówczas otrzymamy po przekształceniu definicji (7): 


Re + |F)7 o” 
v) ny 


Łatwo stwierdzić, że wyrażenie pod pierwiastkiem jest nowym modułem 
bezwymiarowym, który można nazwać zmodyfikowaną liczbą Reynoldsa dla 
mieszania 


Re, e. (8) 


Uwzględniając definicję (8) w równaniu (5) otrzymamy po przekształceniu 
Nu "ly Re; PP" PS: (9) 


Równanie (9) ściślej uzależnia wnikanie ciepła od hydrodynamiki procesu 
aniżeli równanie (5). Stała C, jest w tym przypadku miarą efektywności 
mieszania dla przypadku wnikania ciepła w mieszalniku. 

Równania (5), (6) oraz (9) przybliżają zadowalająco dane doświadczalne 
wnikania ciepła w stosunkowo wąskim przedziale zmian liczb Reynoldsa 1 
Prandtla. Dla bardzo szerokich zmian tych modułów — co ma miejsce 
szczególnie wtedy, kiedy mamy do czynienia z dużą zmianą skali aparatu 
(parametry D, d) oraz zmianą lepkości cieczy — lepszą dokładność uzyskuje 
się stosując do opisu danych doświadczalnych wnikania ciepła równania typu 


3.2 
I 
le [05 | 


3.0 d i 7 QQ 9 
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Nu 
p, 03 V 044 
DN 
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Ryc. 4. Zależność lg(Nu/Pr**> V":'*) =f(Re): 1) mieszadło turbinowe tarczowe, 2) mieszadło 
turbinowe o łopatkach prostych a = 90”, 3) mieszadło śmigłowe 
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Ryc. 5. Wartości poprawek Y, w równaniu (6): a) dla mieszadła turbinowego tarczowego, b) dla 


mieszadeł: turbinowego o łopatkach pochylonych oraz śmigłowego (i — inwariant podobieństwa 
geometrycznego: i, — inwariant dla geometrii standardowej) 
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wielomianów potęgowych. Równania takie nie są jednak wygodne w zapisie 
oraz ponadto nie mają wyraźnego sensu fizycznego. 

Badania doświadczalne wnikania ciepła w mieszalnikach wykonane przez 
różnych autorów [4—6] pozwoliły na wyznaczenie stałych C oraz Cp, jak 
również poprawek Y, dla różnych mieszalników. 

Dobra zgodność wyników istnieje w odniesieniu do wykładników 4, B 
oraz E. Dla cieczy niutonowskich oraz burzliwego przepływu cieczy w 
mieszalniku przyjmuje się [4, 5, 9, 13]: A=0,67; B=0,33 oraz E=0.,14. 

Na rycinie 4 przedstawiono przykładowo wyniki badań doświadczalnych 
współczynnika wnikania ciepła dla trzech typów mieszadeł szybkoobroto- 
wych, zbiorników różnej wielkości oraz burzliwego zakresu przepływu cieczy 
niutonowskiej, uogólnione w postaci funkcji (5). 

Najlepszą zgodność zależności (5) z wynikami doświadczeń uzyskuje 
się [2] dla A =0,67. 

Wartości stałych C oraz C, równań (5) i (6) przedstawiono w tabeli 1, 
gdzie podane są również parametry geometryczne układu mieszadło — 
zbiornik, stosowane w badaniach. Z tabeli tej widać znaczne różnice w 
intensywności wnikania ciepła osiągane przy użyciu różnych mieszadeł. Naj- 
niższą wartość C uzyskano dla mieszadła śmigłowego, najwyższą zaś dla 
mieszadła turbinowego tarczowego o średnicy d = 0,5D oraz dużej liczbie 
łopatek Z = 10. 

Wartości poprawek *, uwzględniających wpływ ważniejszych parametrów 
geometrycznych układu na intensywność wnikania ciepła pokazano na ry- 
cinach 5a oraz 5b. Z rycin tych widać stosunkowo silny wpływ takich 
parametrów jak: wysokość zbiornika (H), liczba łopatek mieszadła (Z) oraz 
wysokość zawieszenia mieszadła (h), natomiast słaby wpływ liczby i szero- 
kości przegród. 

Bardziej wyraźny obraz efektywności pracy różnych mieszadeł dla wnika- 
nia ciepła przedstawia wykres funkcji (9). Na rycinie 6 przedstawiono wyniki 
pomiarów wnikania ciepła dla różnych mieszalników, w tym również zaopa- 
trzonych w niepełne przegrody. Najwyższą efektywność wykazuje mieszalnik 
z mieszadłem turbinowym tarczowym o średnicy d=0,5D i przegrodami 
niepełnej długości [11], przy czym liczba przegród o szerokości B =0,17D 
wynosiła J =8 oraz liczba łopatek mieszadła Z = 10. Najniższą efektywność 
działania wykazuje mieszadło Śmigłowe z przegrodami standardowymi. 


WNIOSKI 


Przeprowadzone badania doświadczalne wnikania ciepła w mieszalnikach 
ujawniły silny wpływ niektórych parametrów geometrycznych układu mie- 
szadło — zbiornik na stałą € w równaniu (5). 

Wykorzystując uzyskane pomiary współczynników wnikania ciepła wy- 
znaczono szereg poprawek Y, występujących w równaniu (6), uogólnionych 
matematycznie, przydatnych do celów projektowych. 


c 


Analiza zależności (9) pozwoliła na stwierdzenie, że różne mieszadła 
wykazują różną efektywność dla przypadku wnikania ciepła w mieszalniku. 

Wykazano, że dla tego przypadku mieszalnik z przegrodami niepełnej 
długości oraz mieszadłem turbinowym tarczowym pracuje najbardziej efek- 
tywnie. 


WYKAZ OZNACZEŃ 


A,B, E[—] — wykładniki w równaniu (5) 

a [m] — długość łopatki mieszadła 

B [m] — szerokość przegrody 

b [m] — szerokość łopatki mieszadła 

C [-] — stała w równaniu (5) 

CG, [-] — stała w równaniu (6) 

Ez L-] — stała w równaniu (9) 

G [J:kg''-K"'] — ciepło właściwe cieczy mieszanej 

D [m] — średnica wewnętrzna mieszalnika 

d [m] — średnica mieszadła 

F [m] — powierzchnia wymiany ciepła 

G [kg's""] — natężenie przepływu cieczy mieszanej 

H [m] — wysokość napełnienia cieczą mieszalnika 

h [m] — odległość dolnej krawędzi mieszadła od dna 
zbiornika 

J L-] — |liczba przegród 

K [-] — stała w równaniu (4) 

k [W:'m'*:K]  — współczynnik przenikania ciepła 

L [m] — długość przegrody 

n [s] -— częstość obrotów mieszadła 

P [W] — moc mieszania 

p [m] — odległość przegrody od dna zbiornika 

Q [W] — strumień ciepła 

q [W :m 7] — jednostkowy strumień ciepła 

Ś [m] — skok mieszadła 

S [m] — grubość ścianki zbiornika 

ty [C] — temperatura cieczy na wlocie do mieszalnika 

t» [C] — temperatura cieczy na wylocie z mieszalnika 

„ad [C] — średnia temperatura cieczy mieszanej 

tp [C] — temperatura pary grzejnej 

| [m] — objętość cieczy mieszanej 

Z [-] — |liczba łopatek mieszadła 

a [”] — kąt pochylenia łopatki mieszadła 

dl [W:-m *:K"'] — współczynnik wnikania ciepła po stronie mie- 


szanej cieczy 


Op [W:m"*:K""] 


7 [kg'm"*'s"*] 
UĄ [kg'm"*:s"'] 
| [W:m ":K"7] 
A, [W:m ":K"'] 
0 [kg'm *] 

rg [4 

; E;1 
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współczynnik wnikania ciepła po stronie kon- 
densującej pary 

dynamiczny współczynnik lepkości cieczy mie- 
szanej 

dynamiczny współczynnik lepkości cieczy mie- 
szanej w temperaturze Ścianki 

współczynnik przewodzenia ciepła cieczy mie- 
szanej 


— współczynnik przewodzenia ciepła ścianki zbior- 


nika 

gęstość cieczy mieszanej 

poprawka zredukowana 

poprawka zredukowana dla i-tego parametru 
geometrycznego 


MODUŁY BEZWYMIAROWE 


D 
Nu = — — |liczba Nusselta 
CH . 
= 7 — |liczba Prandtla 
d 
Re = A : > — liczba Reynoldsa dla mieszania 
(PIV) D*o* , : : 
Re, = > — zmodyfikowana liczba Reynoldsa dla mieszania 
Ą 
U : tyż 
J= k, — simpleks lepkości 


FRYDERYK STRĘK 


EXPERIMENTS ON HEAT TRANSFER IN AGITATED VESSELS 


Summary 


In the article research results of heat transfer in agitated vessels with varying types of 
high-speed stirrers are presented. The measuring was carried out in a turbulent range of the 
liquid's flow in a tank. In the experiments various Ntwtonian fluids and four tanks with fiat 
bottom and diameters varying from 0,15 m to 0,9 m, were used. 

On the basis of simultaneous measurements of heat transfer and intensity of mixing, the 
effectiveness of mixing various mixers was evaluated and the best solution for a given case 


of mixing was suggested. 
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WPBIAJPbBIK CTPEHK 


UCCJIETOBAHHA TEIUIOOTĄNAUM B ANIIAPATAX 
C MEIIIAJJIKAMNM 


Pe3oMe 


B cTaTbe npeACTABJIEHbI pe3yJIbTATbI COOCTBEHHbIX MKCCJIETOBAHHŃ TEIJIOOTTAUH B arnrnapa- 
Tax C pAa3HbIMM THIIAMM ObICTPOBDPALIAFOLNMXCA MELIIAJIOK. FI3MepeHnd BbINOJIHEHbI B OOJIACTH 
TypÓyJIeHTHOTO HBH>KXCHHA HPIOTOHOBCKHX XKMJIKOCTEŃ. VICCJIENOBAHHA IMPOMZBENEHBI B UETbipeX 
cocyH1ax mMaMeTpoM OT 0,15M no 0,9M, € IJIOCKUMH NHMIIAMM H OTpaXKaTEJIbHBIMH Ieperopo1- 
KAMA. 

Ha OCHOBAHHH OJIHOBpeMeHHbIX M3MepeHHi TENIOOTHTAaUM M MOLIHOCTH IIepeM€LUIHBAHHA 
NaHa ONneHka 3dyheKTUBHOCTH paOOTbl pa3HblX ArIapaTOB C MEINAJIKAMM M IIDpENIOXEHbI Halń- 
JyuliMe peLIeHHA MIA TAKOTO CJIYUAA IIEpEMELIMBAHHA. 
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JERZY BOCZAR 


UPROSZCZONE MODELE ROZPRZESTRZENIANIA. 
I PRZEMIAN ZANIECZYSZCZEŃ 
"W ZBIORNIKACH WODNYCH 


Instytut Inżynierii Wodnej Politechniki Szczecińskiej, 
Al. Piastów 50, 70-311 Szczecin 


Streszczenie. Przedstawiono krótkotrwałe zmiany jakości wody, powodo- 
wane przepływami wiatrowymi w zbiornikach wodnych, na przykładzie jeziora 
Miedwie. Dla opracowania prognoz zmian jakości wody na etapie koncepcji i 
założeń techniczno-ekonomicznych zaproponowano stosowanie uproszczonych 
modeli matematycznych rozprzestrzeniania i przemian zanieczyszczeń w zbiorni- 
kach. Modele prędkości przepływu opracowano na podstawie opisu matematy- 
cznego dwuwymiarowego pola prędkości przepływu wiatrowego po uwzględnieniu, 
że znaczący jest gradient prędkości w kierunku prostopadłym do brzegu zbiornika. 
Uproszczony model stężeń zanieczyszczeń opracowano ograniczając go do stężeń 
maksymalnych oraz przez aproksymowanie gaussowskiego rozkładu stężeń wykre- 
sami trójkątów i trapezów. Dla ułatwienia obliczeń przedstawiono nomogramy 
podstawowych wzorów i praktyczne przykłady ich zastosowania. 


Słowa kluczowe: modele matematyczne — rozprzestrzenianie zanieczy- 
szczeń — zbiorniki wodne — przemiany zanieczyszczeń — przepływy wiatrowe — 
pole prędkości przepływu — pole stężeń zanieczyszczeń. 


Jakość wody w zbiornikach określana jest w okresach homotermii głów- 
nie na podstawie Średnich stężeń zanieczyszczeń w całej objętości wody 
zbiornika. W okresach stratyfikacji gęstości, jakość wody opisuje się oddziel- 
nie w warstwie epilimnionu i hypolimnionu. W uzasadnionych przypadkach 
przeprowadzany jest podział zbiornika na wyodrębnione obszary, różniące 
się budową morfologiczną, warunkami hydraulicznymi i sposobem użytko- 
wania wody. 

Średnie stężenia zanieczyszczeń w całej objętości zbiornika lub w jego 
wyodrębnionych częściach określane są na podstawie punktowych pomiarów 
stężeń zanieczyszczeń, dokonywanych najczęściej w miesięcznych przedzia- 
łach czasowych. Na podstawie wyników pomiarów opracowuje się charakte- 
rystykę zmian jakości wody w cyklu rocznym. 
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Przy uwzględnieniu zmian średnich stężeń zanieczyszczeń opracowuje się 
długoterminowe prognozy jakości wody w zbiorniku jako wynik bilansu 
dopływów i odpływów wody, ładunków zanieczyszczeń oraz biochemicznych 
procesów przekształcania zanieczyszczeń w zbiorniku. Długoterminowe prog- 
nozy jakości odnoszą się do cyklu rocznego albo do wieloletnich horyzontów 
czasowych. W „prognozach tych uwzględnia się między innymi takie czynniki 
kształtujące jakość wody, jak elementy eutrofizacji, osady denne, ilość 1 
jakość wodnych organizmów roślinnych i zwierzęcych oraz projektowane 
ładunki zanieczyszczeń odprowadzanych do zbiornika. 

Praktyczne przykłady wykorzystania wody zbiorników do celów gospo- 
darczych, a zwłaszcza poboru wody dla wodociągów komunalnych i potrzeb 
przemysłu, oraz równocześnie dla odprowadzenia wód zanieczyszczonych 
wykazują potrzebę uwzględnienia (w opisie stanu istniejącego i prognozach) 
krótkoterminowych zmian jakości wody. Przyczyną istotnych i częstych 
zmian jakości wody w miejscu jej poboru do celów gospodarczych jest 
dynamika wód zbiornika, w tym głównie przepływy wiatrowe, tranzytowe 1 
cyrkulacyjne oraz — przy większych wymiarach zbiornika — sejsze, genero- 
wane różnicą ciśnienia atmosferycznego. 

Wymienione przepływy powodują rozprzestrzenianie zanieczyszczeń od 
miejsc dopływów wód z ładunkami zanieczyszczeń do miejsc poboru wody, 
zmieniając w sposób znaczący jej jakość w krótkich przedziałach czasowych. 
Zmiany jakości mają istotny wpływ na sposób uzdatniania wody. W przy- 
padkach nieuwzględnienia zmian są one przyczyną zakłóceń w eksploatacji 
urządzeń uzdatniających i w użytkowaniu wody dla gospodarki komunalnej i 
przemysłu. 

Istotne znaczenie w zmianach jakości wody mają dryfowe przepływy 
wiatrowe. Przy prędkościach przepływu, dochodzących do kilkudziesięciu 
centymetrów na sekundę, czas konwekcji zanieczyszczeń nie przekracza 
jednej doby nawet przy kilkunastokilometrowych odległościach pomiędzy 
dopływem zanieczyszczeń i miejscem poboru wody. Przy gospodarczym 
wykorzystaniu zbiornika przez wielu użytkowników czasy przepływu zanie- 
czyszczeń w strumieniach wiatrowych są na tyle krótkie, że procesy bioche- 
miczne nie odgrywają praktycznego znaczenia w zmniejszeniu stężeń zanie- 
czyszczeń. Dwukierunkowe przepływy dryfowo-gradientowe [5] powodują 
rozprzestrzenianie biernych zanieczyszczeń z warstw powierzchniowych do 
przydennych, i odwrotnie, zmniejszając skuteczność selektywnych ujęć wody. 

Przepływy tranzytowe przez zbiornik znaczących dopływów rzecznych i 
przepływy cyrkulacyjne występujące w wyniku ich oddziaływania [6] powo- 
dują wprawdzie rozprzestrzenianie zanieczyszczeń ze znacznie mniejszymi 
prędkościami od wiatrowych, ale w czasie działania wiatru mogą przyczynić 
się do skrócenia czasu przepływu zanieczyszczeń. W strumieniach tranzyto- 
wych występuje znaczący ładunek zanieczyszczeń dopływających ze zlewni. 
Mniejsze znaczenie praktyczne mają przepływy wód zbiornika wywołane 
zjawiskiem sejsz jeziornych. 


49 


Potrzebę korzystania z krótkoterminowych prognoz zmian jakości wody 
w użytkowaniu gospodarczym zbiornika potwierdziły między innymi warun- 
ki eksploatacji ujęcia i stacji uzdatniania wody z jeziora Miedwie, zaopatru- 
jących Szczecin w wodę pitną. 

W dalszej części opracowania przedstawiono praktyczne sposoby obliczeń 
rozprzestrzeniania zanieczyszczeń w zbiornikach wodnych, przy uwzględnie- 
niu dryfowych przepływów wiatrowych wody. Przedstawione sposoby okreś- 
lenia prędkości przepływów wiatrowych i zmiany stężeń zanieczyszczeń sta- 
nowią uproszczenie wcześniejszych, bardziej szczegółowych opracowań anali- 
tycznych autora [3, 5]. Porównanie wyników szczegółowych modeli matema- 
tycznych z wynikami uproszczonych modeli pozwoliło na oszacowanie wiel- 
kości błędu i przedstawienie wystarczająco dokładnych modeli dla praktyki 
inżynierskiej. 


CHARAKTERYSTYKA I ZMIANY JAKOŚCI WODY 
JEZIORA MIEDWIE 


Dla przedstawienia istotnego oddziaływania dynamiki wód, a w tym 
głównie przepływów wiatrowych na dużą częstotliwość zmian jakości wody 
w zbiorniku posłużono się wynikami pomiarów niektórych wskaźników 
jakości wody jeziora Miedwie (ryc. 1, 2). Pomiary przeprowadzane były, na 
ogół codziennie, w trakcie sześcioletniej eksploatacji ujęcia i stacji uzdatnia- 
nia wody dla wodociągu komunalnego w Szczecinie, o wydajności istniejącej 
120 tys. [m*/d] i około 200 tys. [m*/d] po rozbudowie. Pobór wody dla 
wodociągu odbywa się z głębokości +16 m poniżej zwierciadła wody. 
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Ryc. 1. Ogólna liczba organizmów w wodzie pobieranej z ujęcia i w wodzie uzdatnionej w roku 
1982 
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Ryc. 2. Stężenia azotu organicznego w wodzie pobieranej z ujęcia w latach 1980-1981 


Jezioro Miedwie, piąte co do wielkości w Polsce, posiada zlewnię o 
wielkości 1033 km*. Przy średnim rocznym stanie zwierciadła wody 14,03 m 
n.p.m., maksymalna długość zbiornika wynosi 16,2 km, szerokość 3,6 km, 
głębokość 44 m (średnia 19,3 m), powierzchnia zwierciadła wody 3527 ha 
oraz objętość 681 672 tys. m”. Rynna zbiornika ma kierunek południkowy. 
Oś podłużna zbiornika pokrywa się z linią łączącą maksymalne głębokości. 
Głównym dopływem jeziora jest rzeka Płonia o średnim rocznym przepływie 
około 4,5 m*/s. Pozostałe dopływy: Ostrawica, Gowienica i Miedwianka 
posiadają łączny przepływ Średni roczny około 0,95 m*/s. 

W okresie letnim występuje okresowa stratyfikacja wód jeziora, z termo- 
kliną na głębokości około 10 do 12 m. Wyraźny brak stratyfikacji występuje 
w okresach jesiennym i wiosennym. Jezioro zamarza całkowicie przy wyjąt- 
kowo mroźnych zimach. Pokrywa lodowa utrzymuje się na ogół w strefach 
przybrzeżnych. 

Jezioro jest eksponowane na działanie wiatrów, o zmieniających się z 
dużą częstotliwością kierunkach i prędkościach od 2 do 15 m/s. Najczęściej 
występuje wiatr z kierunku SW z prędkością 5 m/s. Istotne dla zmian jakości 
wody w miejscu ujęcia dla wodociągu są działania wiatrów z kierunków SE 1 
NE, powodujące dopływy wody z rejonu ujścia rzeki Płoni i akwenu 
kąpieliska i ośrodka rekreacyjnego w Morzyczynie. 

Poza poborem wody ze zbiornika dla zaopatrzenia Szczecina w wodę 
pitną jest on wykorzystany do celów rekreacyjnych oraz przez gospodarstwa 
rybackie. Z głównego kąpieliska i obiektów rekreacyjnych korzysta maksy- 
malnie około 20 tys. osób dziennie. Małe obiekty wypoczynkowe i żeglarskie 
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usytuowane są również na zachodnim brzegu zbiornika. Istotnymi źródłami 
zanieczyszczenia wody zbiornika są dopływy rzek. 

Do rzeki Płoni odprowadzane są Ścieki komunalne i przemysłowe z 
Barlinka, Pyrzyc i Kłodzina (przy czym tylko w Barlinku znajduje się 
sprawna oczyszczalnia Ścieków) oraz z pięciu gorzelni rolniczych. Ze względu 
na intensywne wykorzystanie zlewni rzeki Płoni i pozostałych rzek do celów 
rolniczych, wody powierzchniowe doprowadzają do jeziora duże ładunki 
substancji biogennych. 

Jak wynika z warunków eksploatacji stacji uzdatniania wody dla wodo- 
ciągu, istotnymi wskaźnikami jakości wody jest ogólna liczba organizmów 
oraz stężenie azotu organicznego (ryc. 1, 2). Dominującą grupę w planktonie 
stanowią Ssinice oraz w okresie wiosennym okrzemki. Ogólna liczba organiz- 
mów w okresie intensywnego rozwoju fitoplanktonu wzrasta do ponad 15 
tys/cm* i stanowi istotne utrudnienie w eksploatacji filtrów pospiesznych 
i uzyskaniu właściwych efektów uzdatniania wody. Stężenia azotu orga- 
nicznego przekraczają bardzo często wielkości dopuszczalne dla pierwszej 
klasy jakości wody powierzchniowej. Przewidywany dalszy rozwój gospo- 
darki rolnej oraz słabe tempo budowy oczyszczalni Ścieków w zlewni jeziora 
wskazują na możliwość zwiększenia ładunku związków azotu i fosforu od- 
prowadzanych do jeziora i intensyfikację procesów eutrofizacyjnych. 

Sporadyczne pomiary wykonane w rejonie czynnego ośrodka rekreacyjne- 
go i ujścia rzeczki Miedwianki wykazały zanieczyszczenie bakteryjne wody. 

Przepływy wiatrowe z kierunku ujścia rzeki Płoni, gdzie występują duże 
stężenia substancji odżywczych i intensywny wzrost fitoplanktonu przenoszą 
zanieczyszczenia w rejon ujęcia wody, zwiększając ogólną liczbę organizmów 
i stężenia azotu organicznego w wodzie. Podobny wpływ mają intensywne 
przepływy wiatrowe dryfowo-gradientowe, które powodują podniesienie osa- 
dów z dna zbiornika. Przepływy wiatrowe z kierunku ośrodka rekreacyjnego 
w Morzyczynie mogą powodować zmiany jakości wody zanieczyszczeniami 
bakteryjnymi. 

Eksploatacja stacji uzdatniania wody dla wodociągu wykazała konie- 
czność prognozy zmian jakości wody, przy uwzględnieniu dynamiki wód 
zbiornika, a w szczególności przepływów wiatrowych. Wielkość zmian i ich 
częstotliwość mają praktyczne znaczenie dla przyjęcia prawidłowej technolo- 
gii uzdatniania wody oraz technicznego wyposażenia stacji. 


UPROSZCZONY SPOSÓB OBLICZANIA ROZPRZESTRZENIANIA 
I PRZEMIAN ZANIECZYSZCZEŃ 


Jakość wody w jeziorach i zbiornikach opisywana jest często przy zasto- 
sowaniu dynamicznych modeli matematycznych. Modele te charakteryzują 
się złożoną strukturą i wymogiem przeprowadzenia pracochłonnych czynnoś- 
ci dla estymacji albo identyfikacji parametrów i weryfikacji modeli. Do tej 
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grupy modeli należą niektóre modele ekologiczne, których rozwój uzasadnio- 
ny jest intensyfikującymi się procesami eutrofizacji zbiorników. 

Ogólna postać modelu ekologicznego opisującego na przykład rozprze- 
strzenianie i przemiany fitoplanktonu w wodzie zbiornika ma postać nastę- 


pującą: 


0S Ta A ME 08 ( 08 
Ot 0x 0y Gz  ćxl *6x] cyl ” óy 


ż ( + x,y,z)'S-Dft, x,y,z):S (1) 


gdzie: 

S — jest stężeniem fitoplanktonu, 

G — jest funkcją opisującą czynniki wzrostu masy fitoplanktonu 
(natężenia Światła w warstwie eufotycznej, temperatury wody, 
stężenia związków biogennych); 

funkcjaD — opisuje czynniki zmniejszania masy fitoplanktonu (obumieranie 
komórek, respirację, ubytek fitoplanktonu spożywanego przez 
zooplankton). 
Wpływ natężenia światła na funkcję wzrostu fitoplanktonu opisany został 
przez Jórgensena [9]: 


j lę(t)+K, 


e jo REL 2 
kzę  I(z,t1) ŁK, (2) 


G[Io(t)] = 


gdzie: 

Ią(t) [kcal:m *:d'"] — jest natężeniem światła na powierzchni zwierciadła 
wody w chwili £t, 

k [m'':s] — jesc wspołczynnikiem pochłaniania światła, 

K;, — jest stałą Michaelisa. 


I(z, t) = Io(t)exp(=kz); (3) 
k=«a+$58 
a [m '] — współczynnik pochłaniania Światła przez wodę; 


B [m*:g''] — współczynnik pochłaniania światła spowodowany zawiesiną 
fitoplanktonu. 

Wpływ stężenia substancji biogennych na funkcję wzrostu fitoplanktonu, 
zgodnie z teorią Monoda, jest następujący: 


b, 
G(b,) = —— (4) 
(b, + b;) t 
gdzie: 
b, — jest stężeniem i-tej substancji biogennej, 
b, — jest stężeniem, przy którym funkcja wzrostu fitoplanktonu jest równa 


połowie maksymalnej wartości funkcji wzrostu w optymalnych warunkach. 
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Wpływ temperatury T na funkcję wzrostu określony jest na podstawie 
zależności [1, 11]: 


G(T) = Kr:T (5) 


gdzie: Kr[d ':*C''] jest stałym parametrem. 
Zmniejszenie masy fitoplanktonu na skutek respiracji zależy od tempera- 
tury i opisane jest jako [2]: 


D(T) = R'h(T) (6) 
gdzie: R[d '] — maksymalna wielkość ubytku fitoplanktonu na skutek 
respiracji, 
oraz 
DJ ao Y:(Tmax 7 Topt) 
h(T) e = "exp [y i (Tys EJ kę, JAR (7) 
JAR e b i 
gdzie: 
LA temperatura, przy której występuje maksymalna wartość funkcji 


wzrostu, 
1,,, — temperatura, powyżej której nie występuje dalsze zwiększenie funkcji 
wzrostu G, 
y — stały parametr. 
Ilość spożytego fitoplanktonu określona jest jako [9]: 


l 
BEŻ EZ (8) 
gdzie: 
Z [g'm *”] — masa zooplanktonu, 
n, [d''] — współczynnik wzrostu zooplanktonu, 
a [g'm *] — współczynnik asymilacji pochłoniętego fitoplanktonu. 


Wielkość funkcji D (zmniejszenia fitoplanktonu na skutek obumierania 
komórek) przyjmowana jest jako wielkość stała [9] bądź zależna od tem- 
peratury. 

Jak wynika z przedstawionych zależności (2--8), ta część modelu (1), która 
dotyczy zmiany masy fitoplanktonu, w wyniku procesów wzrostu i zmniejsze- 
nia wymaga estymacji i identyfikacji parametrów na drodze pracochłonnych 
pomiarów. 

Stąd równanie (1) w ogólnej postaci, uwzględniającej przemiany i rozprze- 
strzenianie fitoplanktonu w wodzie zbiorników, stosowane jest rzadko [11]. 
W większości stosowanych modęli uproszczono ponadto opis transportu 
masy fitoplanktonu, ograniczając się do opisu wymiany pionowej pomiędzy 
warstwami epilimnionu i hypolimnionu, lub w przypadku dużych zbiorników 
do opisu transportu poziomego pomiędzy wydzielonymi obszarami zbior- 
nika, przyjętymi jako jednorodne [2, I, 9, 13]. 
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W większości modeli albo nie uwzględnia się rozprzestrzeniania masy fito- 
planktonu w wyniku przepływów wiatrowych, albo przyjmuje się pośrednie 
oddziaływanie wiatru na zmniejszenie stratyfikacji wody w zbiorniku [8]. 
Dotyczy to okresów stratyfikacji wody, kiedy w powierzchniowej warstwie 
może występować wysoka temperatura wody, sprzyjająca rozwojowi fito- 
planktonu. Przy turbulencji wody wywołanej działaniem wiatru następuje 
częściowa redukcja stratyfikacji i obniżenie temperatury wody, a tym samym 
zmniejszenie masy fitoplanktonu. Zjawisko odwrotne wzrostu masy fito- 
planktonu występuje podczas intensywnej cyrkulacji wiatrowej 1 falowania, 
które mogą spowodować okresową likwidację termokliny w zbiorniku 
1 wzrost ilości związków biogennych zgromadzonych w hypolimnionie, rów- 
nież w warstwie osadów. 

Poniżej przedstawia się opis matematyczny rozprzestrzeniania, który 
uwzględnia przepływ wody wywołany działaniem wiatru. 

Przy strumieniu wiatrowym, o kierunku zbliżonym do wyrównanej linii 
brzegowej zbiornika (ryc. 3f), rozprzestrzenianie zanieczyszczeń odbywa się 
przy oddziaływaniu powierzchni zwierciadła wody, dna i brzegu zbiornika. 
Oddziaływanie brzegu zbiornika można pominąć, jeżeli strumień wiatrowy 
oddala się od brzegu (ryc. 3g). 

W opisie matematycznym uwzględniono brzegowe wyloty ścieków, usy- 
tuowane: 1) przy powierzchni zwierciadła wody (ryc. 3a—c): bezpośrednio 
przy brzegu (ryc. 3a) albo w odległości „b” od brzegu (ryc. 3c), albo 2) w 
odległości „h* od powierzchni dna (ryc. 3d—e). Przedstawione przypadki 
usytuowania wylotu ścieków odpowiadają najczęściej występującym w prak- 
tyce. 

Pole stężeń zanieczyszczeń przy założeniu ustalonego przepływu wiatro- 
wego może być opisane równaniem konwekcji-dyfuzji [3] 


08 08 08 08 


= [854 (a, Bai Ł.S 5 
Toe SAO AEP Ay SAME BÓL SE ©) 

Po określeniu pola wektorów prędkości przepływu wiatrowego na pod- 
stawie układu równań (28 -- 30), po estymacji współczynników dyfuzji burzli- 
wej 1 współczynnika szybkości przemian zanieczyszczeń niestabilnych [4], 


przy założeniu ustalonego dopływu ścieków i przy uwzględnieniu warunków 
brzegowych 


sę (10) 


układ równań (9, 10) może być aproksymowany równaniami roznicowymi 
przy niejawnym schemacie różnicowym i rozwiązany numerycznie [3]. 
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: WE ES ABOCA V=VZM=ZME 


Ryc. 3. Schematy rozprzestrzeniania zanieczyszczeń: f) wzdłuż brzegu, g) od brzegu zbiornika, 
i usytuowania wylotu Ścieków (a --e) 


Że względu na konieczność zagęszczenia dyskretnych punktów siatki 
obejmującej pole stężeń zanieczyszczeń, obliczenia przeprowadzane są przy 
zastosowaniu e.m.c. 

Poniżej przedstawia się propozycję uproszczonego sposobu określenia 
rozprzestrzeniania zanieczyszczeń, który może być wykonany przy zastoso- 
waniu prostych obliczeń 1 pomocniczych wykresów. Dokładność wyników 
obliczeń jest wystarczająca do wstępnej prognozy oddziaływania odprowa- 
dzanych zanieczyszczeń na jakość wody w zbiorniku, zwłaszcza w miejscach 
poboru wody do celów gospodarczych. 

Zakłada się przepływ jednowymiarowy w strumieniu wiatrowym z jedna- 
kową prędkością v, w kierunku równoległym do kierunku działania wiatru. 
Przestrzenne pole stężeń zanieczyszczeń może być określone na podstawie 
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analitycznego rozwiązania równania (9) po uwzględnieniu przepływu jedno- 
wymiarowego (v, = v, =0) oraz po założeniu dopływu zanieczyszczeń, usy- 
tuowanego według schematu „a” na ryc. 3 [3]: 


B B 
21,——-y acz 


S(x,Y,Z)= "e Jet (—— ) erf | x 
> 


X 
2 [Bt Sa 
H 
577 2,——+z 
<|et( ser — | |araa (11) 
ZD 2 Dy 


gdzie: 
wyrazy w nawiasie kwadratowym — reprezentują konwekcję 1 dyfuzję burzli- 
wą zanieczyszczeń, 


współczynnik a, — oddziaływanie brzegu, 

współczynnik «a, — oddziaływanie powierzchni zwierciadła wody i dna 
zbiornika, 

współczynmk a, — zmniejszenie stężenia zanieczyszczeń niestabilnych w 


wyniku procesów przekształcania. 
Uproszczone obliczenia ogranicza się do maksymalnego stężenia w prze- 
kroju poprzecznym strumienia. Przy współrzędnych 


l 
sza [srt(—F |rert(- : ) |--a. (12) 
NA Y) X ; X 


p. Dy = 
v [2 


y 


Ta sama zależność jest ważna przy usytuowaniu wylotu ścieków przedsta- 
wionym na schematach „b=e” na ryc. 3. 
Oznaczając konwekcję i dyfuzję zanieczyszczeń w przestrzeni nie ograni- 


czonej symbolami (5; 1 (Z) , zamiast zależności (12) otrzymuje się 
0/y 0/z 


= (5 | (5 


W dalszych obliczeniach rozprzestrzeniania uwzględniono zanieczyszcze- 
nia bierne, co oznacza, że gęstość i iepkość ścieków nie różni się od wody. 
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Tym samym współczynnik burzliwej dyfuzji zanieczyszczeń w kierunku pio- 
nowym D, jest równy współczynnikowi burzliwej wymiany pędu. 


; X ; > ; ; 
Przy czasie przepływu — zanieczyszczeń na drodze od wlotu do zbiornika 
U 


do miejsca poboru wody, nie przekraczającym kilkunastu godzin, wynik 
procesu przemian zanieczyszczeń jest na tyle mały, że pozwala na przyjęcie 
wielkości a, = 1. 

Wymiary powierzchni początkowej przy punktowych wylotach ścieków 
oblicza się na podstawie zależności 


źli el 2 = Li (14) 


Przy dopływach rzecznych wymiary 2/, i 21, odpowiadają średniej szerokości 
i głębokości koryta rzeki na odcinku ujściowym. Dla ułatwienia obliczeń 
według wzoru (13) opracowane zostały wykresy (ryc. 4). 

Wielkości współczynników a, i «, mogą być określone na podstawie 
warunku ciągłości ładunku zanieczyszczeń. Przy założeniu dwuwymiarowego 
procesu rozprzestrzeniania w nieograniczonej płaszczyźnie x, y równanie 
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ciągłości ma postać [4] 
ON ODC (15) 
0 
Powierzchnia nieograniczonęgo pola stężeń zanieczyszczeń w dowolnym 


przekroju, przy zaniedbaniu procesów przemian, wynosi [4] 


je ©) 


s 1 I-y © 
M5 [|er( SD |ror(-* ] || = 255 l,. (16) 
o 2 |D,— 2 |D,Ę 


a D 


Powierzchnia pola stężeń ograniczonego brzegami zbiornika 


rl W 


B 
Dz) 


a, 
Ś "2 AA 
A> = — erf Rap = +erf = =F dy+ 
2 Xx X 
= 2 |D,> 2 |D,> 
50 | B B 
MSEZARW ZAL | 
Sp 2 2 
+ 7 erf RE W =—==—galii x xi erf SE z dy : ( 17) 
X X 
ZE, 2 |D,= 2 |P,; 


Współczynnik uwzględniający oddziaływanie brzegów a, może być określony 
z zależności (16, 17) 


A ud 4l, 
y A, z 
s =Y zTY 
| er = rat( = |o+ 
Ę 5 1pŹ ep 
PG *v "e 
(18) 
_B 
2 2 = PRA 
F | jex( — rat = 
L X x 
Bl, 2 D,> 2 D,> 


Podobnie jak (18) może być określony współczynnik a,, uwzględniający 
oddziaływanie powierzchni dna 1 zwierciadła wody. Całki w mianowniku (18) 
obliczane są przybliżoną metodą całkowania, przy zastosowaniu programu 
na maszynę cyfrową [4]. 

Dla uproszczonych obliczeń wielkości współczynników a, i «,, uwzględniają- 
cych oddziaływanie powierzchni zwierciadła wody, dna i brzegów zbiornika 
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proponuje się, aby pola stężeń opisane krzywą Gaussa (16, 17) zastąpić 
powierzchnią trójkątów i trapezów. Jak wynika z przeprowadzonego porów- 
nania, wielkości współczynników a, i a, obliczone metodą uproszczoną 
różnią się na tyle nieznacznie od wielkości obliczonych według (18), że mogą 
być stosowane w obliczeniach praktycznych maksymalnych stężeń. 

Współczynnik a, dla wylotu usytuowanego bezpośrednio przy brzegu 
(ryc. 3a) może być określony w sposób uproszczony z porównania nieograni- 
czonego pola stężeń w dowolnym przekroju z polem stężeń ograniczonym 
brzegiem zbiornika (ryc. 5a). 


Ryc. 5. Schematy pola stężeń zanieczyszczeń 


Powierzchnia sprowadzona nieograniczonego pola stężeń zanieczyszczeń, 
w przekroju poprzecznym, gdzie znany jest stosunek stężenia maksymalnego 


do stężenia w wylocie |— | ścieków wynosi 
0/y 


a|S 
A, =-|—| =4dl. 
: (5) ; (19) 


Powierzchnia pola stężeń zanieczyszczeń z uwzględnieniem oddziaływania 


brzegu zbiornika 
ar |2rfa l 
mek >) > 2 
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a przy uwzględnieniu (19): 


r SED 24 
2 (5) AB 8 (21) 
So/y 
Z warunku ciągłości A; = A> otrzymuje się: 
r 4 Ś 
WE Ej 13 a (5, G 
y 


So), 8150 
Współczynnik a, dla tego samego wylotu ścieków (ryc. 3a) określić można 
podobnie przy uwzględnieniu schematu „b” (na ryc. 5): 


S S S 
rH (>) rH 3 rl, (5) 
i So So z So z (23) 


y 


81, 4 4 


H_ 4 
kt = : ZE 1 w (ÓŃ) 


Sa s) s WE ie H* Ly kt 
So: Sos So). 4. 8I3 4 Śo/ż 


Po podstawieniu do (24) H = th otrzymuje się zależność (22). 


A, =rH— 


Ta sama wielkość (24) współczynnika a, występuje w przypadku wylotu 
ścieków oddalonego od brzegu na odległość „b” (ryc. 3c) oraz przy rozprze- 
strzenianiu zanieczyszczeń w strumieniu przepływu wiatrowego, oddalającego 
się od brzegu (ryc. 3b). 

Współczynnik a, dla wylotu usytuowanego w odległości „b” od brzegu 
(ryc. 3c—e) określono na podstawie schematu „d” (ryc. 5): 


4 b. 4 
FR Zs, MŚ 


e E (+? Ę EA) Pi 
PSE E AŚCYGA8" 41, =8I7 zw 


dla b = O otrzymuje się zależność (22). 
Wielkość współczynnika a«,, przy wylocie ścieków usytuowanym w odleg- 
łości „h” nad dnem zbiornika (ryc. 3d —e), określono na podstawie schematu 
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Po podstawieniu do (26): H AA h=0 


otrzymuje się zależność (22). 
Zależności (22) oraz (24—26) i ich wykresy (ryc. 6; 7) stosuje się przy 
obliczeniu stosunku maksymalnych stężeń zanieczyszczeń przez pomnożenie 


[S 
wielkości (5) i (5) „ określonych dla nieograniczonej przestrzeni 
0/y 0/z 


przez odpowiedni współczynnik a, i a,. 
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Ryc. 6. Współczynnik oddziaływania brzegów zbiornika a, według wzoru (25) 


Maksymalne wielkości iloczynów spełniają ponadto warunek 


Ś 
Ę ByPSSLO: 
S0 y,z | 


Wielkość współczynnika zmniejszającego stężenie niestabilnych zanieczy- 
szczeń a, określa się na podstawie znanej zależności 


4, = exp (-k 3 (27) 


Współczynnik szybkości zmiany stężenia zanieczyszczeń k,[s '] w wy- 
niku ich przekształcenia może być okresiony na podstawie pomiarów tere- 
nowych. W uproszczonych obliczeniach maksymalnych stężeń zanieczyszczeń 
(13) w strumieniach wiatrowych, w przypadku kiedy czas przepływu od 
wylotu ścieków do miejsca poboru wody wynosi od kilku do kilkunastu 
godzin, można przyjąć odpowiednio 


OK = LO = 0,9. 
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Ryc. 7. Współczynnik oddziaływania dna i zwierciadła wody a, 


Pola wektorowe prędkości przepływów wiatrowych określa się na podsta- 
wie znanego układu równań ruchu burzliwego cieczy nieściśliwej oraz równa- 
nia ciągłości. Pomijając w równaniach ruchu wpływ przyspieszenia Coriolisa 
i pionowe składowe prędkości, które mają znikome znaczenie dla prakty- 
cznych obliczeń przepływów wiatrowych w zbiornikach wodnych, otrzymuje 
się 

OÓv, Óv, Óv, lóp 10, 


FSA NOŚ af 00 po AJK x 2 
o 0x ”óy  0dx o Óz S9 


— HV ŁU 5 ———— — (29) 
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Ga: b) 0; (30) 

Óx  0y 
Rozwiązanie układu równań (28 -30) przeprowadza się przy założeniu usta- 
lonych przepływów wiatrowych, bez uwzględnienia stratyfikacji gęstości wo- 
dy na całej głębokości zbiornika w okresach homotermii (wiosna, jesień), 
albo w warstwie epilimnionu w pozostałych okresach roku. Po wprowadze- 
niu do układu (28 —30) funkcji wydatków jednostkowych Y(x, y) otrzymuje 
się równanie różniczkowe eliptyczne, które wraz z warunkami brzegowymi 


F|r = 0; m p, = 0; P| rs m 
gdzie: 
I, — stały brzeg zbiornika, 
I, — wodna granica zbiornika, 
Ia — dopływ wód rzecznych do zbiornika, 


można rozwiązać numerycznie, przy zastosowaniu niejawnego schematu róż- 
nicowego. 

W opracowaniach [3, 5] przedstawione są rozwiązania tego równania dla 
stałego i zmiennego współczynnika dyfuzji, wraz z przykładami obliczeń. 

Wymienione modele szczegółowe charakteryzują się złożonością struktu- 
ry, pracochłonnością i kosztownością środków obliczeniowych. Na ich pod- 
stawie, po sprawdzeniu wielkości błędu, mogą być opracowane modele 
uproszczone, korzystniejsze dla stosowania w praktyce inżynierskiej biur 
projektowych i administracji gospodarki wodnej. 

Poniżej przedstawia się uproszczony sposób obliczenia Średniej prędkości 
jednowymiarowego przepływu wiatrowego, który może być zastosowany w 
praktyce przy wykorzystaniu prostych wzorów i wykresów, bez potrzeby 
stosowania programu na maszynę cyfrową. 

Po uwzględnieniu w równaniach (28, 30), że przepływ wiatrowy jest 
ustalony 


że w strefie przybrzeżnej zbiornika znaczące zmiany pola prędkości prze- 
pływu występują w kierunku prostopadłym do brzegu [7] (ryc. 8) 


A asa" 


oraz że: 
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Ryc. 8. Układ współrzędnych do opisu pola prędkości: I — wyrównane dno zbiornika, 2 — 
statyczne zwierciadło wody, 3 — położenie zwierciadła wody podczas przepływu wiatrowego 


układ równań (28 --30) przyjmuje postać: 


v, 06 10, 
o = Maz eż S1) 
Óv, _ OŚ 100, 
ia td, wózo 02 


Analityczne rozwiązanie układu równań (28 —30) [12], przy założeniu 
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Po podstawieniu do równań (33, 34) 
TOx 2 i Toy kn ; ZAS: ; nn -—6 
—- =qHo' Wiocos0; —-=4o Wyosin6; 49=3:10", 
0 0 
stałego współczynnika dyfuzji burzliwej pędu 
D.=q4-WięH. gdzie: «4 =0,54=10* (35) 
i współczynnika oporów ruchu przy dnie 
Na = 8,5-102-n2-H73 (36) 


gdzie: n — jest współczynnikiem Manninga, 
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oraz po założeniu, że kierunek wiatru jest zgodny z kierunkiem osi układu 
(ryc. 8), 
dla wiatrów wiejących od brzegu w kierunku zbiornika albo w kierunku 
odwrotnym otrzymuje się: 

1 


W, Ź 2 || W 2 z? 
py = o (+4 +35 - (2 (1735 | (37) 
da H* 


HH 


dla wiatrów równoległych do wyrównanej linii brzegu: 
1 


W z Wh 
z (1 był | i (38) 


NH 


Pionowy rozkład prędkości w przepływie dryfowym równoległym do brzegu 
(ryc. 9a) uwzględnia wpływ tarcia na powierzchni zwierciadła wody i dna 
zbiornika. 


w W 


Ryc. 9. Rozkład prędkości w pionie: a) przepływ wzdłuż brzegu, b) przepływ od brzegu 


W strumieniu płynącym w kierunku od brzegu zbiornika występuje, poza 
przepływem dryfowym, gradientowy przepływ — w wyniku gradientu ciśnie- 
nia hydrostatycznego, wywołanego obniżeniem zwierciadła wody przy brzegu 
(ryc. 9b). W miarę oddalania się od brzegu przepływ gradientowy maleje, a 
dryfowy wzrasta, w wyniku czego pionowy rozkład prędkości przybliża się 
do rozkładu przedstawionego na rycinie 9a. Kolejna zmiana rozkładu pręd- 
kości występuje na przeciwległym brzegu zbiornika. 

Dla określenia pola stężeń zanieczyszczeń przyjęto rozkład prędkości jak 
w strumieniu dryfowym, który występuje podczas przepływu wzdłuż brzegu 
oraz przeważa w przypadku przepływu w kierunku od brzegu zbiornika. 

Prędkość średnia w strurnieniu dryfowym wynosi (ryc. 10): 


c H 


5 — Szczecińskie Rocz. Nauk. t. IL. z. I 
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Tabela 1. Współczynnik n, dla obliczenia prędkości wiatru na wysokości 10 m nad powierzchnią 
zwierciadła wody 


Charakterystyka terenu w miejscu | ż. 


Lp. Nye M AW LIWA 
położenia stacji pomiarów prędkości wiatru B* >5km | B<5 km | 


| |W lesie, wiatrowskaz na poziomie wierzchołków drzew, 
pojedyncze wierzchołki powyżej wiatrowskazu 3,0 12 


2 |W lesie, wiatrowskaz znacznie powyżej wierzchołków 
drzew e) 1,6 


3 | Na granicach zabudowy miejskiej, pojedyncze budynki albo 

drzewa w odległości 30 --40 m (z jednej strony) przewyższają | 

wiatrowskaz | 2,0 1,3 | 
4 | Jak w p. 3, wiatrowskaz znacznie powyżej drzew i budynków | 1,8 1,1 | 


5 |W osiedlu wiejskim, na odkrytej, równej płaszczyźnie; 
najbliższa zabudowa w odległości 30--50 m, wiatrowskaz 


powyżej drzew i budynków 1,7 1,05 
6 | Na odkrytej płaszczyźnie (pole, łąka, lotniska), 

iw odległości 200 --500 m las albo zagajnik 1,6 1,0 
7 | Na brzegu jeziora, morza albo większej rzeki, w odległości 

100 --200 m las albo budynki 1,4 0,9 
8 | Na równym terenie niezadrzewionym 1,3 0,8 
9 | Na całkowicie odkrytym brzegu jeziora, morza albo 


większej rzeki, bezpośrednio przy zwierciadle wody ji! 0,7 


10 | Na całkowicie odkrytej, wąskiej mierzei, bezpośrednio 
przy zwierciadle wody 1,05 0,65 


* szerokość zbiornika 


Tabela 2. Współczynnik n, dla obliczenia prędkości wiatru na wysokości 10 m nad powierzchnią 
zwierciadła wody 


Stromy wierzchołek wzgórza, przy urwisku 
Zaokrąglony wierzchołek wzgórza, górna krawędź zbocza 


Równina, bardzo szeroka dolina 


Dolna część zbocza, dno wąskich i płytkich dolin, kotlin, wąwozów 


Dno głębokich dolin, kotlin, wąwozów 
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Wp [ms"*] 
100 —. — + DZ a | 
80—— R= | ŻA | 
6,0 | Z | | 
s_- Ą 
MAE 
RSMA v [ems*] 


Ryc. 10. Średnie prędkości przepływu wiatrowego = Wio(zg+ jąc) 1) Hy = 0,05; 
UJ! 


2) Nu 5 0,10; 3) ny = 0,20; 4) ny = 0,40; 5) ny = 0,55. 


6 18 20 


Ryc. 1l. Współczynnik oporów ruchu przy dnie zbiornika qy =8,5:10*:n*:H"*3% 


Współczynnik oporów ruchu przy dnie może być określony dla współczynni- 
ków Manmninga: 

n =Q0,016 — dno gładkie piaszczyste, 

n=Q0,018 — zwarty żwir na dnie, 

n =0,020 — dno ilaste, 
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n =Q0,025 — dno pokryte roślinnością (na podstawie wykresu — ryc. 11). 
Przy zmieniającej się powierzchni dna na odcinkach I, drogi przepływu 
wiatrowego można zastosować zastępczy współczynnik tarcia przydennego 


> Nail ż 
HH,z — = ] | . (40) 


Prędkość wiatru na wysokości 10 m nad powierzchnią wodną zbiornika 
określa się na podstawie prędkości pomierzonej na stacji meteorologicznej, 
pomnożonej przez współczynniki korygujące [10] n;, i n, (tab. 1, 2). Współ- 
czynniki uwzględniają charakterystykę terenu (n,) oraz rzeźbę terenu (n,) 
w miejscu położenia stacji meteorologicznej. Współczynnik (w) uwzględnia 
wpływ prędkości wiatru na iloczyn współczynników korygujących: 


URZYERUE 


STĘŻENIA ZANIECZYSZCZEŃ W JEZIORZE MIEDWIE 


Obliczenia stężeń zanieczyszczeń wody przeprowadzono dla przypadku 
przepływów wiatrowych wzdłuż zachodniego brzegu jeziora. 

W pierwszym przykładzie obliczono maksymalne stężenie ogólnej liczby 
organizmów transportowanych z rejonu ujścia rzek Płoni i Ostrowicy do 
miejsca ujęcia wody dla wodociągu, na drodze o długości około 11,0 km. 
W południowej .części jeziora występują sprzyjające warunki do rozwoju 
fitoplanktonu, w wyniku znacznego ładunku substancji biogennych, dopły- 
wających głównie ze zlewni rzeki Płoni. 

Prędkość wiatru z kierunku S i SE przyjęto: W =8,0 m/s na podstawie 
danych z najbliższej stacji meteorologicznej usytuowanej w Szczecinie-Dąbiu, 
w odległości około 30 km. 

Tę samą prędkość wiatru przyjęto do obliczeń przepływu wiatrowego po 
uwzględnieniu charakterystyki i rzeźby terenu w miejscu położenia stacji 
pomiarów oraz prędkości wiatru: n, = 1,0; n, =1,0;: w=L0 (tab. 1—3). 
Średnią głębokość wody na drodze przepływu wiatrowego: H =4,5 m, 
przyjęto na podstawie wymiarów ośmiu przekrojów poprzecznych strefy 
przybrzeżnej o szerokości około 30 m. 

Współczynnik oporów ruchu przy dnie uwzględnia warunki dna pokrytego 
roślinnością (n = 0,025). Na podstawie zależności (36) i wykresu (ryc. 11) 


1 
Hy = 8,5107 -0,0252-4,5 3 =0,32. 


Prędkość średnią w strumieniu dryfowym otrzymuje się na podstawie zale- 
żności (39) i wykresu (ryc. 10): v = 0,248 [m's"']. 
Współczynnik burzliwej dyfuzji zanieczyszczeń biernych w kierunku piono- 
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Tabela 3. Współczynnik w, uwzględniający prędkość wiatru 


wym określono na podstawie zależności (35) 
D. =0,54* 10 *:80-45= 1.97410% [ms 2]. 


Estymacja współczynnika dyfuzji zanieczyszczeń w kierunku poziomym, pro- 
stopadle do głównego kierunku rucnu, przeprowadzana była na podstawie 
pomiarów terenowych stężeń rodaminy w sześciu przekrojach poprzecznych 
strumienia. Wielkość współczynnika obliczona metodą Fischera wynosi: 


Du 16877" m 2 gl 


Powierzchnię początkowego przekroju poprzecznego określono przy 
uwzględnieniu całej szerokości strefy przybrzeżnej 2/, = 15 m oraz intensyw- 
nego rozwoju organizmów w warstwie powierzchniowej o grubości równej 
zda =lg" m 

Stosunek maksymalnych stężeń S ogólnej liczby organizmów do stężeń 
początkowych S$, obliczono na podstawie zależności (12, 13) 1 wykresu 
(ryc. 4), 


x 11000 
la - = = 44 500 
"p o 0,248 ; 
7,5 
(5; = erf (—— | =e(1393 
DY 2.4407 3-44500, 
(3) = er zm 2) 550) = 0,029 
So): 2 ,/1,97:107*-44500 


po wyznaczeniu «, i w, na podstawie zależności (22, 24) i wykresów (ryc. 6, 7) 
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otrzymano 


a, = 1,575; a, = 15,9. 


Przyjmując, że w stosunkowo krótkim czasie przepływu strumienia wia- 
trowego — = 12,4 godz. występuje zbliżone tempo rozwoju i ubytku ogólnej 
v 


liczby organizmów w wyniku procesów biochemicznych (a, = 1), otrzymuje 
się w miejscu ujęcia wody: 


S 
ię 0,573 : 1,575-0,029:15,9 = 0,416. 


0 


Jeżeli zatem w rejonie ujścia rzeki Płoni, przy sprzyjających warunkach, 
występują maksymalne ilości organizmów o stężeniu S, = 30 tys. [cm *], 
to w wyniku przepływów wiatrowych mogą być one przeniesione w rejon 
ujęcia wody dla wodociągu, przy zmniejszeniu stężenia w wyniku dyfuzji 
burzliwej — do około 12,5 tys. [cm '*]. 

W drugim przykładzie obliczeń określono maksymalne stężenia zanieczy- 

szczeń dopływających w rejon ujęcia wody dla wodociągu podczas przepływu 
wiatrowego od ujścia rzeczki Miedwianki, usytuowanego w rejonie ośrodka 
rekreacyjnego w Morzyczynie. Odległość wzdłuż brzegu zachodniego wynosi 
około 5,8 km. 
Prędkość występujących wiatrów z kierunków N i NE określono jak po- 
przednio: W =8 [m:s"']. Średnia głębokość wody na drodze przepływu 
wiatrowego obliczona w czterech przekrojach poprzecznych strefy przybrze- 
Żnej o szerokości 20 m wynosi 5,2 m. Współczynnik oporów ruchu przy dnie 
wynosi na podstawie (36) i wykresu (ryc. 11) 


Ny = 4,5-102-0,0252-5,27 13 = 0,30. 


Prędkość średnią przepływu wiatrowego określono na podstawie (39) 
i wykresu (ryc. 10) v = 0,25 [ms "]. 

Współczynnik dyfuzji burzliwej zanieczyszczeń w kierunku pionowym okreś- 
lono na podstawie (34): 


D, = 0,54-107*-8,0-5,2 = 2,28-1073 [m?2-s"1] 


oraz w kierunku poziomym na podstawie pomiarów terenowych 


D,= 1073 [m?:s"*]. 


Powierzchnia przekroju początkowego wypływu zanieczyszczonych wód 
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rzeki Miedwianki wynosi 2l, = 3,0 m; 21, =0,5 m 


REES 


b 023 


1,50 
(5; R er >) — 0,174 
SA 24/1073-23200 


0,25 
(5) = er „ya )- 0018. 
Sg: 2, 228-10 *-23200 


Oddziaływanie brzegu oraz powierzchni dna i zwierciadła wody na proces 
rozprzestrzeniania uwzględniają współczynniki a, = 1,84; a, = 11,30, określo- 
ne na podstawie (22, 24) 1 wykresów (ryc. 6, 7). Przy zaniedbaniu zmniejsze- 
nia stężeń w wyniku procesów biochemicznych (a, = 1), co uzasadnia stosun- 


kowo krótki czas przepływu Ę = 6,45 godz.) otrzymuje się stosunek maksy- 
malnego stężenia w miejscu ujęcia wody (S$) do stężenia początkowego w 
miejscu wypływu zanieczyszczeń do jeziora: 

a = 0,174: 1,84 :0,018: 11:30 = 0,065. 

So 

Mając na względzie wymagane w miejscu ujęcia wodociągowego dopu- 

szczalne stężenie zanieczyszczeń wody klasy pierwszej: 

BZT. 44 [mg'dcm *], a miano coli > 1, 

dopuszczalne stężenia zanieczyszczeń, w miejscu ich wylotu do jeziora winny 
mieścić się w normie: 

BZT; £< 61 [mg'dcm*”] oraz miano coli > 0,065. 

Mając na względzie miano coli, jakość wody w rejonie ośrodka rekreacyj- 
nego i kąpieliska powinna być lepsza, niż to określają przepisy dla drugiej 
klasy jakości, wymaganej do urządzenia zorganizowanych kąpielisk (miano 
coli 2 0.1). 


WNIOSKI 


Przepływy wody w zbiorniku wywołane działaniem wiatru mają istotny 
wpływ na zmiany jakości wody. Powodują bowiem rozprzestrzenianie zanie- 
czyszczeń dopływających do zbiornika z wodami rzek i przez wyloty Ścieków, 
doprowadzając je do miejsc ujęć wody dla celów komunalnych i przemysło- 
wych. Potwierdza to analizowany przykład jeziora Miedwie. Wynika z tego, 
że źródła zanieczyszczeń usytuowane w odległościach około 11 km (dopływ 
rzeki Płoni z dużym ładunkiem substancji biogennych) 1 około 6 km (rejon 
ośrodka wypoczynkowego w Morzyczynie i dopływ rzeki Miedwianki) powo- 
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dują na skutek przepływów wiatrowych znaczące i częste zmiany jakości 
wody w miejscu ujęcia dla wodociągu miejskiego. Czas przepływów wiatro- 
wych, przy wymienionych odległościach, jest na tyle krótki, że zmniejszenie 
stężeń zanieczyszczeń w procesach biochemicznych nie ma praktycznego 
znaczenia. 

Dla niezawodności działania urządzeń ujmujących i uzdatniających wodę 
zbiorników wymagane jest opracowanie poza długoterminowymi prognozami 
zmian jakości wody w okresie roku lub wielolecia także krótkoterminowych 
prognoz zmian jakości wody w wyniku przepływów wiatrowych. 

Dla opracowania krótkoterminowych prognoz zmian jakości wody mogą 
być zastosowane uproszczone modele rozprzestrzeniania 1 przemian zanieczy- 
szczeń. Uproszczone modele opracowane zostały na podstawie wykonanych 
wcześniej szczegółowych modeli matematycznych pól prędkości przepływów 
wiatrowych i pól stężeń zanieczyszczeń w zbiornikach wodnych. 

Uproszczony analityczny sposób obliczeń pola prędkości przepływów 
wiatrowych otrzymuje się przez uwzględnienie w układzie równań dwuwy- 
miarowego ruchu burzliwego i równania ciągłości przepływu, rozwiązywa- 
nych numerycznie. Pozwala on stwierdzić. że w strefie przybrzeżnej zbiornika 
znaczące zmiany prędkości przepływu 1 różnice poziomów zwierciadła wody 
występują w kierunkach prostopadłych do brzegu. 

Uproszczony sposób obliczenia dwuwymiarowego pola stężenia zanieczy- 
szczeń polega na ograniczeniu algorytmu do maksymalnych stężeń w prze- 
krojach poprzecznych strumienia oraz aproksymowania, przyjętego w szczegó- 
łowych modelach rozkładu zanieczyszczeń według krzywej Gausa, rozkładem 
trójkątnym albo trapezowym. 

Poza uproszczeniem opisu matematycznego i algorytmu, w przedstawio- 
nych modelach zastosowano proste wzory analityczne i opracowane na ich 
podstawie nomogramy (zamiast analitycznych albo numerycznych rozwiązań 
za pomocą programów na maszyny cyfrowe). 

W przedstawionych przykładach zastosowania uproszczonych modeli wy- 
niki obliczeń oceny zmian jakości wody są wystarczająco dokładne na etapie 
opracowań koncepcyjnych i założeń techniczno-ekonomicznych ujęć 1 stacji 
uzdatniania wody ze zbiorników. Prognozy krótkoterminowych zmian jakoś- 
ci wody generowanych głównie przepływami wiatrowymi, opracowane wed- 
ług uproszczonych modeli, mogą stanowić razem z prognozami zmian długo- 
terminowych podstawę określenia warunków niezawodnego działania urzą- 
dzeń do ujmowania i uzdatniania wody zbiorników. 


WYKAZ OZNACZEŃ 


— szerokość zbiornika [m] 
— odległość wylotu zanieczyszczeń od brzegu zbiornika [m] 
xyz  - współczynnik dyfuzji burzliwej [m*:s'"] 


OFTU 
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g — przyspieszenie ziemskie [m:s 7] 

H — głębokość wody [m] 

h — odległość wlotu zanieczyszczeń od dna zbiornika [m] 
k, — współczynnik szybkości przemian zanieczyszczeń [s '] 
ży szerokość powierzchni początkowej [m] 

21, — wysokość powierzchni początkowej [m] 

l, — długość i-tego odcinka dna zbiornika [m] 

n — współczynnik oporów ruchu Manninga 

p — ciśnienie hydrostatyczne [N'm ?] 

Q, — dopływ rzeczny [m*'s "] 

Q, — dopływ wody zanieczyszczonej [m*'s "] 

S — stężenie zanieczyszczeń [Gm *]; [mg'dm *] 

NY — stężenie zanieczyszczeń w wylocie ścieków [G:m *] 

I — czas [s] 

Dyys  — prędkość przepływu wody [m's '] 

Wio — prędkość wiatru na wysokości 10 m nad powierzchnią zwierciad- 


ła wody [m's '] 

x,y,z — współrzędne układu [m] 

pr — granica lądowa zbiornika 

Ć — wzniesienie dynamicznego zwierciadła wody ponad poziomem 
statycznym [m] 


Hp — współczynnik tarcia przy dnie zbiornika 

UŃ — współczynnik aerodynamicznego tarcia na powierzchni zwierciad- 
ła wody 

o — gęstość wody [kg'm ”] 

aj — naprężenia styczne na powierzchni zwierciadła wody [N'm 7] 

dk: — naprężenie styczne [N'm 7] 


JERZY BOCZAR 


SIMPLIFIED MODELS OF POLLUTION PROPAGATION 
AND TRANSFORMATION IN WATER RESERVOIRS 


Summary 


Short — term changes in the quality of water — resulting from wind flows in water 
reservoirs are presented on the example of Miedwie Lake. To work out prognoses for changes in 
the quality of water at the stage of conception, technical and economical assumptions, it is 
suggested that simplified mathematical models of pollution propagation and transformation in 
reservoirs be used. Models of flow velocity have been worked out on the basis of mathematical 
description for a two-dimensional velocity field of wind flow. It had been taken under account 
that the velocity gradient perpendicular to the bank of the reservoir is the significant one. The 
simplified model of pollution concentration has been worked out by limiting it to maximal 
concentrations and by approximating the Gaussian distribution of concentrations with graphs of 
triangles and trapezia. To facilitate calculations, nomograms of basic formulas and examples of 
their application are presented. 
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EXXM BOUAP 


YIIPOIIEHHBIE MONEJIM PACIIPOCTPAHEHH A 
A IPEBPALIEHUM 3ATPA3HEHAMHM B BOHIOEMAX 


Pe3oMme 


IIpencTaBJIeHbl KpaTKOBpeMEHHbiE M3MEHeHKA KAUECTBA BOJĄbI, BbIZBAHbIe BETDOBbBIMH TeUe- 
HHMAMM B BOJTOĆMAX. Ha rpumMepe oj3epa MeBe. /Jl1a pazpaó0Tku nporHo30B HU3MEeHEHNA 
KAUECTBA BOJIbI Ha 3TAIMe KOHUEMIMM HM TeXHMHECKO-JKOHOMMUECKOTO IIOJJIOKEHMA IpeJUIo- 
>XCHO IIDpAMEHCHHE YIDPOLIGHHbIX MATeMATHUECKHX MOJeJleń pacHpocTpaHeHH4 M rnpeBpauieHuń 
3aTpA3HeHHA B BOJOĆMAX. MOJNEJIM CKODOCTH IIpoTekaHHA OOÓpaGOTaHbl Ha OCHOBe MATEMaTH- 
4eCKOTO OIIACAHHA TOJIA CKODOCTM BETPpOBOTo TeUeHHA, KOTOpPOE HUMEET NBA M3MEPEHHA, C 
yHUGTOM TOTO, UTO 3HAUMMbIM ABJIAETCA TpANMEHT CKOPOCTM, HarNpaBJIeHHbIA repreH MKyJlAaDHO 
k 6epery BO1OĆMa. 

Y rpolieHHyio MOJIEJIb KOHIIEHTDANNMM 3ATpAZHCHMHA BbIIIOJIHMJIM OTDAHMUKBAA EĆ JO MAK- 
CHMAJIBHbIX KOHIIEHTpaliMi, a TAaKYKXe METOHOM AIIHpoOKCHMAIIMA pacnpenejleHuA KOHIEHTPALMM 
no Iayccy, maarpamMMaMn TpeyTOJIbHMKOB M Tparelnń. JlJIa ynpolieHna pacHćTOB Ipe1CTaB- 
JIEHbl MOHOTDAMMbi OCHOBHbIX (DODMyJI M IpaKTMUECKKE IpMMEpbl MX IpHMEHCHMA. 
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Streszczenie. W pracy przedstawiono metodę kształtowania właściwości 
dynamicznych obrabiarek. Oparta jest ona na współzależności zagadnień modelo- 
wania, identyfikacji i modyfikacji. Identyfikację parametrów modeli można prze- 
prowadzić na podstawie niepełnej informacji z badań doświadczalnych. Do minima- 
lizacji kryterium identyfikacji stosowane są metody gradientowe. Algorytm modyfi- 
kacji umożliwia wielokryterialną ocenę rozwiązań. Jako zmienne decyzyjne w 
procesie modyfikacji przyjmuje się te parametry konstrukcyjne, które mają naj- 
większy wpływ na przyjęte kryteria oceny właściwości dynamicznych. Optymalne 
wartości zmiennych decyzyjnych wyznacza się metodą systematycznego przeszuki- 
wania punktów planu eksperymentu numerycznego. Zaprezentowana metoda two- 
rzenia planów eksperymentu umożliwia równomierne rozmieszczenie zadanej licz- 
by punktów w obszarze ograniczeń. Przedstawiono przykład praktycznego wyko- 
rzystania opracowanego algorytmu kształtowania właściwości dynamicznych. 


Słowa kluczowe: obrabiarka — dynamika — modelowanie — identyfika- 
cja — modyfikacja. 


WPROWADZENIE 


Prognozy rozwojowe obrabiarek skrawających do metali jednoznacznie 
wskazują na wzrost znaczenia procesów dynamicznych zachodzących w 
układzie OUPN (obrabiarka — uchwyt — przedmiot obrabiany — narzędzie 
skrawające) [1]. Wynika to stąd, że obrabiarki najbliższej przyszłości charak- 
teryzować się będą znacznie większymi prędkościami obrotowymi i posuwo- 
wymi oraz znacznie większymi mocami silników napędowych. Jest to deter- 
minowane przede wszystkim rozwojem w dziedzinie materiałów narzędzio- 
wych, a także ciągłym dążeniem do podnoszenia wydajności przy równoczes- 
nym wzroście dokładności, trwałości użytkowej i niezawodności działania. 
Rozpatrywana w takich uwarunkowaniach jakość użytkowa obrabiarek bę- 
dzie w znacznym stopniu zależna od właściwości dynamicznych układu 
OUPN. Problemem o pierwszorzędnym znaczeniu staje się wytwarzanie 
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obrabiarek o właściwościach dynamicznych spełniających wymagania wynm- 
kające z przedstawionych wyżej trendów rozwojowych. 

Kształtowanie właściwości dynamicznych w procesie projektowania obra- 
biarki nie może być dziełem przypadku, lecz powinno być działaniem w pełni 
świadomym. Konstruktor powinien zatem być „wyposażony” w odpowiedni 
algorytm kształtowania właściwości dynamicznych obrabiarki. Obecny stan 
wiedzy z dziedziny dynamiki obrabiarek umożliwia opracowanie efektywnego 
algorytmu, który pozwoli na praktyczną realizację zadania kształtowania 
właściwości dynamicznych układu OUPN już na etapie projektowania obra- 
biarki. 

W procesie tworzenia nowej konstrukcji dokonuje się odwzorowania 
zbioru parametrów (geometrycznych, materiałowych itp.) opisujących kon- 
strukcję w zbiór wskaźników służących do oceny jakości użytkowej obra- 
brarki (ryc. 1) 


A 
LP) > U) (1) 
gdzie: 
lp — zbiór parametrów opisujących konstrukcję, 
(U! — zbiór wskaźników oceniających jakość użytkową obrabiarki, 
A  — operator przyporządkowujący elementom zbioru |p) elementy zbio- 
ru (U. 


Zbiór wskażników 
jakości opracowanej 
konstrukcji (u) 
Podzbiór 
wskaźników 
jakości dyna- 
micznej 


Zbiór struktur 

i parametrów 
opisujących 
konstrukcję | pył 


Ryc. 1. Odwzorowania w procesie kształtowania właściwości dynamicznych obrabiarek 


Odwzorowanie to ma charakter zwężający, bardzo liczny bowiem zbiór 
parametrów |p| opisujących konstrukcję zostaje odwzorowany w mało liczny 
zbiór (U) wskaźników jakości. Odwzorowanie to może odbywać się w 
sposób świadomy (konstruktor w procesie konstruowania korzysta z dostęp- 
nych metod obliczeniowych, pozwalających na wyznaczenie wartości liczbo- 
wych wskaźników jakości) lub przypadkowy (konstruktor w procesie kon- 
struowania kieruje się intuicją, wzoruje się na rozwiązaniach istniejących). 
Uzyskany przez odwzorowanie proste (wzór 1) wynik może znacznie odbie- 
gać od wyniku o pożądanych właściwościach. Korzystając z odwzorowania 
prostego konstruktor może rozwiązania pożądanego poszukiwać tzw. metodą 
prób i błędów. Metoda ta jednak nawet przy znacznej liczbie iteracji nie 
gwarantuje uzyskania rozwiązania pożądanego, a tym bardziej optymalnego. 

Z. tego względu, z punktu widzenia konstruktora, bardziej racjonalna jest 
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możliwość praktycznego przeprowadzenia odwzorowania odwrotnego zbioru 
wskaźników w zbiór parametrów 


(U) — (pl. (2) 


Zgodnie z wzorem (2) zbiór wskaźników o wartościach odpowiadających 
pożądanym właściwościom obrabiarki (w tym także dynamicznym) jest od- 
wzorowany za pomocą operatorów A ' w zbiór parametrów opisujących 
konstrukcję. Odwzorowanie (2) charakteryzuje się niejednoznacznością prze- 
jawiającą się w tym, że jednemu elementowi zbioru (U! można przypisać 
nieskończenie wiele elementów zbioru |p|. Dodatkową komplikacją przy 
próbach efektywnego przeprowadzania odwzorowania (2) jest nieznajomość 
analitycznej postaci operaiorów A '. W tej sytuacji praktyczne przeprowa- 
dzenie odwzorowania odwrotnego (2) jest możliwe przy użyciu operatorów 
A ' w postaci odpowiednich algorytmów. 


ALGORYTM KSZTAŁTOWANIA WŁAŚCIWOŚCI DYNAMICZNYCH 
UKŁADU OUPN 


W Instytucie Technologii Mechanicznej Politechniki Szczecińskiej opraco- 
wano algorytm umożliwiający efektywne wykonanie odwzorowania odwrot- 
nego przy rozwiązywaniu zagadnienia kształtowania właściwości dynami- 
cznych obrabiarki [2, 9]. Schemat blokowy tego algorytmu przedstawiono 
na ryc. 2. Idea algorytmu oparta jest na współzależności zagadnień modelo- 
wania, identyfikacji 1 modyfikacji (optymalizacji). Podstawą algorytmu jest 
model obliczeniowy (matematyczny), budowany w oparciu o rysunek kon- 
strukcyjny (złożeniowy), odzwierciedlający z określoną dokładnością właści- 
wości dynamiczne obrabiarki (lub dowolnego jej zespołu). Parametry struktu- 
ralne p, (bezwładnościowe, dysypacyjne, sztywnościowe) są funkcjami para- 
metrów konstrukcyjnych p, (geometrycznych, materiałowych itp.) 


Ps = U (px) - (3) 


W praktyce obliczeniowej najczęściej stosowane są modele liniowe, których 
typowymi przedstawicielami są: 
w przypadku tłumienia wiskotycznego 


Mj+ Hąq+ Kq =P, (4) 
w przypadku tłumienia konstrukcyjnego 
l 
Mj+—Lq+Kq=P, (5) 
©) 
gdzie: 


M — macierz współczynników bezwładności, 
H — macierz współczynników tłumienia wiskotycznego, 
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Ryc. 2. Schemat blokowy algorytmu kształtowania właściwości dynamicznych obrabiarek 


L — macierz współczynników tłumienia konstrukcyjnego, 

K — macierz współczynników sztywności, 

ij, 4, q — wektory (macierze kolumnowe) przyśpieszeń, prędkości i przemie- 
szczeń uogólnionych, 

P — wektor (macierz kolumnowa) uogólnionych sił wymuszających. 


W wyniku rozwiązania modelu matematycznego uzyskuje się charakterystyki 
dynamiczne q (częstotliwościowe, czasowe), które są funkcjami parametrów 
strukturalnych modelu i odpowiednio częstotliwości lub czasu 


di (0) ra! Ja (Dx» w) (6) 
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lub 
di (£) WE Ja (P;, t). (7) 


Na podstawie charakterystyk (6) i (7) wyznacza się wartości wskaźników U 
przyjętych do oceny. 


U =J[4:(0)] = 4 Ufa LM (Px)» o]; = F, (PĄ, 0) (8) 
lub 
U = f5[q,(1)] = fs Us LA (P), t]) = Fa(Py, 1). (9) 


Wzory (8) 1 (9) przedstawiają odwzorowanie proste. Postać analityczna 
operatorów F;, F, jest praktycznie znana tylko w przypadku prostych (o 
małej liczbie stopni swobody) modeli matematycznych. W przypadku modeli 
złożonych operatory F;, F, stosowane są w praktyce obliczeniowej w postaci 
odpowiednich algorytmów. 

Obliczone wartości wskaźników U porównuje się z wartościami dopu- 
szczalnymi lub pożądanymi. Dopuszczalne lub pożądane wartości wskaźni- 
ków mogą być ustalone na podstawie norm państwowych, branżowych, 
zakładowych, na podstawie wartości wskaźników wyznaczonych dla obra- 
biarki uznanej jako wzorcowa lub mogą być ustalone arbitralnie przez 
konstruktora. Jeśli porównanie wypadnie korzystnie dla opracowywanej 
konstrukcji i konstruktor nie chce poszukiwać lepszego rozwiązania — 
można na tym obliczenia zakończyć. Jeśli porównanie wypadnie niekorzyst- 
nie lub konstruktor zdecyduje się na poszukiwanie rozwiązania lepszego, 
wówczas przystępuje do poszukiwania parametrów, które decydują o ocenia- 
nych właściwościach dynamicznych rozpatrywanej konstrukcji. Parametry te 
wyznaczyć można przy użyciu metody wrażliwości. Metoda ta pozwala 
określić, które spośród wszystkich parametrów opisujących model mają 
największy wpływ na wskaźniki oceny. Zmianę AU wskaźnika jakości w 
funkcji zmian Ap; parametrów opisujących model wyznaczyć można korzy- 
stając z zależności: 


Og; (0) 
AU = AD: (9 
A Ops; 
lub 
Oq, (t) 
a DE aa (11) 
gdzie: 
dg. Og, (t 
m, We — to funkcje wrażliwości pierwszego rzędu. 
Dsi Psi 


Znak funkcji wrażliwości wskazuje kierunek, w jakim należy dokonać zmiany 
wartości parametru p,. Jako zmienne decyzyjne wybiera się zatem te para- 
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metry p,,, dla których funkcje wrażliwości przyjmują największe wartości. 
Wybrane zmienne decyzyjne podlegają analizie pod kątem praktycznej reali- 
zowalności ich zmian, po czym ustala się ostatecznie podzbiór |p$; parame- 


trów decyzyjnych oraz określa się dopuszczalne granice ich zmian 


d Z. | 
DSi min 3 Dsi ag? 


spełniających warunek nienaruszenia cech funkcjonalnych obrabiarki. Gra- 
nice wartości tych parametrów wyznaczają obszar, w którym poszukuje się 
rozwiązania pożądanego (optymalnego). 

Przy poszukiwaniu rozwiązania pożądanego (optymalnego) stosuje się 
metodę systematycznego przeszukiwania z planowaniem eksperymentu nu- 
merycznego. Dla obszaru wyznaczonego przez 


d o K: 
Psi min > Psi mx: 


tworzy się plan eksperymentu numerycznego, w którym liczba poziomów 
zmienności każdego parametru jest równa liczbie planowanych doświadczeń 1 
każdemu poziomowi odpowiada jedno doświadczenie. Równomierność roz- 
mieszczenia punktów w przestrzeni parametrów uzyskuje się przez minimali- 
zację kryterium 


l 
DRB2.. 2, a (12) 


gdzie: 
l, — odległość między i-tym a j-tym punktem w unormowanej przestrzeni 
parametrów, 
m — liczba planowanych doświadczeń. 
Na rycinie 3 pokazano plan eksperymentu dla dwóch parametrów decyzyj- 
nych i 20 doświadczeń (m = 20). 

W każdym z punktów planu eksperymentu oblicza się wartości wskaźni- 
ków oceny. Wskaźniki U, podlegają odpowiedniej normalizacji niezbędnej ze 
względu na obiektywizację oceny za pomocą kryterium globalnego Ig. 


Ryc. 3. Plan eksperymentu numerycznego dla 2 zmiennych i 20 doświadczeń 
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Normalizację przeprowadza się w taki sposób, aby wartość każdego ze 
wskaźników oceny zawierała się w przedziale <0, 1). 
Jedynym z możliwych sposobów normalizacji jest nieliniowe odwzorowanie: 


U. U; +e, (13) 
c "gi U, +h, 


Współczynniki przekształcenia (13) wyznacza się z warunków: 
U, =aa>Uxy =0 
U, =beUy=sm 0<m<l (14) 
U =G>Uy=l. 


Na rycinie 4 pokazano przykład normowania wskaźników oceny w 


Ryc. 4. Przykład normowania wskaźników oceny 


przestrzeni dwuwymiarowej. Punkt A — oznacza rozwiązanie idealne, B — 
rozwiązanie oczekiwane, € — rozwiązanie najgorsze, ale jeszcze dopuszczal- 
ne. 

Ponieważ oceny właściwości dynamicznych projektowanej obrabiarki do- 
konuje się na podstawie wielu wskaźników (wielokryterialna ocena rozwią- 
zań), niezbędne jest określenie kryterium globalnego 

ie. Ig(Un,» UN: ++ U)» 
jednoczącego zbiór ocen (U N;j> jakimi scharakteryzowane jest rozwiązanie w 


jedną ocenę całościową. 
Jako kryterium globalne można zastosować jedną z funkcji 
k 
a= ck AU: IG= [l Ux dę= Ż. A; UG; Ją = max U, (15) 
gdzie: A, — odpowiednio dobrane współczynniki masy. 


6 — Szczecińskie Rocz. Nauk. t. II, z. I 
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Wybór postaci kryterium, jak również ustalenie wartości współczynników 
masy nie jest łatwe. Proponuje się ustalenie postaci kryterium globalnego na 
podstawie stopni kompromisu, tj. możliwości kompensacji braku jednych 
właściwości nadmiarem innych [2]. W takim przypadku jako kryterium 
globalne można przyjąć funkcję 

Ee 

le=|ĘŻ 05 |: (16) 

k i=1 | 
Przez odpowiedni dobór parametru r można wtedy uwzględnić różne przy- 
padki zależności ocen od wartości wskaźników. Wartość parametru r określa 
się na podstawie średniego stopnia kompromisu 4 


m Ink 


r = ; (17) 
a. 
In (+ ') 
m 
przy czym 
1 k 
EAC (18) 
ki=ą 
oraz 
Ol; — U, —m, (19) 
gdzie: 
0, — stopień kompromisu dla i-tego wskaźnika, 


Ux, — wartość unormowanego i-tego wskaźnika w punkcie d, (ryc. 4), 
m — ocena rozwiązania w stosunku do ocen 0,1 (0<m<1). 

Za optymalne (w sensie przyjętych wskaźników i kryterium globalnego) 
uważa się te wartości zmiennych decyzyjnych, dla których wartość kryterium 
globalnego jest najmniejsza. Po wyznaczeniu optymalnych wartości parame- 
trów (zmiennych decyzyjnych) sprawdza się, czy uzyskana poprawa właści- 
wości dynamicznych jest zadowalająca oraz czy realizacja optymalnych para- 
metrów jest praktycznie możliwa. Jeśli nie, wówczas należy przeanalizować 
ograniczenia nałożone na parametry decyzyjne, i jeśli jest to uzasadnione, 
dokonać odpowiednich zmian ich wartości. Tworzy się nowy plan ekspery- 
mentu numerycznego, wyznacza się wartości wskaźników, przeprowadza ich 
normalizację i oblicza wartości kryterium globalnego. Parametry odpowiada- 
jące minimalnej wartości kryterium globalnego są parametrami optymalnymi. 
Dokonuje się sprawdzenia, czy poprawa właściwości dynamicznych jest 
zadowalająca i czy realizacja parametrów optymalnych jest praktycznie 
możliwa. Może się okazać, iż mimo wielokrotnego przejścia tej pętli iteracyj- 
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nej nie uzyska się pożądanych właściwości dynamicznych w sensie przyjętego 
kryterium globalnego. W takiej sytuacji można przyjąć najlepsze spośród 
wszystkich rozpatrywanych rozwiązań albo dokonać bardziej radykalnych 
zmian w przyjętym rozwiązaniu konstrukcyjnym (np. w zakresie struktury). 


ALGORYTM IDENTYFIKACJI STRUKTURALNYCH PARAMETRÓW 
MODELI 


Dokładność obliczeń związanych z kształtowaniem właściwości dynami- 
cznych jest ściśle zależna od dokładności, z jaką zastosowany w tych 
obliczeniach model matematyczny będzie odzwierciedlać właściwości dynami- 
czne obiektu rzeczywistego (obrabiarki). Błędy w określeniu wartości parame- 
trów strukturalnych opisujących model (współczynniki bezwładności, tłumie- 
nia, sztywności) oraz błędy w określeniu rodzaju i struktury modelu powodu- 
ją, iż jest mało prawdopodobne, aby zbudowany tylko na drodze obliczeń 
analitycznych model odzwierciedlał z praktycznie wystarczającą dokładnością 
właściwości dynamiczne obiektu rzeczywistego. Zatem jeśli jest to możliwe 
(istnieje np. prototyp obrabiarki), należy zweryfikować model przez porówna- 
nie charakterystyk obliczonych z charakterystykami wyznaczonymi doświad- 
czalnie. Jeśli występują istotne rozbieżności między tymi charakterystykami, 
konieczne jest przeprowadzenie identyfikacji strukturalnych parametrów mo- 
delu. Identyfikację tę można przeprowadzić zgodnie z algorytmem przedsta- 
wionym na ryc. 2. 

Obliczone na podstawie liniowego dyskretnego modelu charakterystyki 
częstotliwościowe 4 = 4(p,, w) są funkcjami parametrów opisujących model. 
Wyznaczone analitycznie parametry strukturalne p, są obarczone błędami. 
Miarą niepewności co do wartości tych parametrów może być macierz 
kowariancji W. Najczęściej macierz ta określana jest intuicyjnie. 

Podstawą do porównania zgodności modelu z obiektem są doświadczal- 
nie wyznaczone charakterystyki częstotliwościowe obiektu — q(w). Charakte- 
rystyki te obarczone są błędami. Błędy systematyczne są korygowane pod- 
czas opracowywania wyników badań [6], natomiast błędy losowe określane 
są przez macierz kowariancji R. 

Że względu na obiektywne trudności pomiarowe oraz koszt badań zwykle 
wyznacza się tylko niektóre charakterystyki i to w ograniczonym paśmie 
częstotliwości. W takiej sytuacji (niekompletna liczba charakterystyk) przyję- 
cie w procedurze identyfikacji jako zmiennych decyzyjnych wszystkich para- 
metrów opisujących model doprowadziłoby do powstania dużych błędów 
obciążenia ich estymat. W opisywanym algorytmie przyjmuje się zatem jako 
zmienne decyzyjne tylko te parametry, które mają największy wpływ na 
kryterium zgodności charakterystyk. 

Do przeprowadzenia oceny wpływu parametrów p,, na charakterystykę 
częstotliwościową q, definiuje się dwa współczynniki wrażliwości [7]: 
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względny współczynnik wrażliwości 


_ 0dilo) _ Odi(0) 


= = Ę 20 
"._ Olnp.j UB mó e 
oraz logarytmiczny współczynnik wrażliwości 
= Olnq, (w) = oq, (w) Psi | (21) 
ln Psj Cpsj d; (w) 
ży: 
Pochodne 4i(0) wyznacza się na podstawie równania 
Psi 
W(o)q(0) = P(0), (22) 


które uzyskuje się w wyniku dokonania transformacji Laplace'a (przy zero- 
wych warunkach początkowych) równań (4) lub (5) i podstawieniu s = jo. 
Różniczkując równanie (22) względem parametru p,, otrzyma się równanie 
macierzowe 


WZA (23) 


Rozwiązując równanie (23) otrzyma się poszukiwane pochodne th 

sj 
Współczynniki I, i I, są funkcjami częstotliwości, zatem wybór parametrów 
do identyfikacji na podstawie analizy ich przebiegu może sprawiać trudności. 
Z tego względu zdefiniowano [8] względne współczynniki wrażliwości pas- 
mowej oceniające globalnie wpływ parametru p,, na charakterystykę q, w 
paśmie częstotliwości 40, Gy). 
Są to: 
— średni względny współczynnik wrażliwości pasmowej 


k-1 
> (0) (On —0,)| 
My m Z (24) 
Ż. 4; (0) (©ną 1 — 0) 


— średni logarytmiczny współczynnik wrażliwości pasmowej 


k— 


1 
MII PY 3 I,(0,) 
n=l 


(0x1 —0,) 
(0, —04) 


(25) 


Parametry p,,, dla których współczynniki MI, i MI, przyjmują największe 
wartości, obiera się jako zmienne decyzyjne ps, w procesie identyfikacji. 
Konieczność wyrażenia zgodności charakterystyk modelu (obliczenio- 
wych) i obiektu (doświadczalnych) w kategoriach liczbowych wymaga przyję- 
cia określonej postaci kryterium identyfikacji. Postać tego kryterium zależna 
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jest od dostępności informacji o niepewnościach zawartych w q oraz p,, od 
przyjętego modelu niepewności oraz od postaci estymatora, który będzie 
stosowany do przekształcenia obserwacji q w estymaty p; rzeczywistych 
wartości parametrów p, . 


Przyjmując bayesowski model niepewności oraz estymator maksimum 
prawdopodobieństwa, wówczas pg wyznaczyć można minimalizując kryte- 
rium 


1 | 
J,(p) => L4—4(pJ] R '[4—4(p)]+[Ps-Psol W  [ps—Psol; (26) 


gdzie: 

p, — wektor parametrów strukturalnych wyznaczonych analitycznie. 
Przyjmując fisherowski model niepewności oraz estymator o największej 

wiarygodności, wówczas pę wyznaczyć można minimalizując kryterium 


|| 


Stosując podejście oparte na minimalizacji ważonej sumy kwadratów nie 
wprowadza się żadnego modelu niepewności. Modelowanie niepewności i 
twierdzenie o optymalności estymatorów zastępuje się intuicyjnym stwierdze- 
niem, że przy danej q sensowną estymatą $ą otrzyma się minimalizując 
kryterium 


1 
J; (p) = 3 Lq—4(pJ]" A14—4(P.)]) (28) 


gdzie: 
A — dodatnio określona macierz wagowa, której elementy wyznacza się na 
podstawie rozeznania inżynierskiego. 

Ponieważ zależność charakterystyk częstotliwościowych q(p,) od parame- 
trów pg jest nieliniowa, zadanie minimalizacji kryterium (we wszystkich 
trzech rozważanych przypadkach) rozwiązuje się przy użyciu iteracyjnych 
metod gradientowych (które cechują się dobrą zbieżnością), jeśli pochodne 
kryterium można wyznaczyć w sposób analityczny. Gradient kryterium wy- 
znacza się na podstawie zależności 


VJ (pw) > — 9 Gl4—4(P.)] (29) 
gdzie: 
b — macierz jakobianowa funkcji wektorowej 4(p,): 
Od, 
Qi, = 0 
Ps; 


G — macierz wagowa, która w zależności od przyjętego modelu niepewności 
jest równa R"" lub A. 
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+ * Proces minimalizacji funkcji J (p,) realizowany jest wg następującej zale- 
Żności rekurencyjnej 


ANA = ps, + Aps, Er Ps, + 44 dy (30) 
gdzie: 
b; — wektor zmiennych decyzyjnych wyznaczonych w k-tej iteracji, 
d,=M VJ (ps) — wektor kierunku poszukiwań, 
V, — macierz modyfikacji kierunku poszukiwań, 


vJ(p.,) — gradient funkcji J(p;) w punkcie Ps,» 
4, — współczynnik kroku minimalizujący funkcję J (Ps) na kierunky poszu- 
kiwań d,. 
Proces iteracyjny prowadzi się do czasu uzyskania zbieżności, to znaczy 
gdy: 
Pb ŻA 
W pracy [2] do minimalizacji kryterium J (p.) stosowano metodę linearyzacji 


modelu Gaussa-Newtona, metodę zmiennej metryki Davidona-Fletchera— 
Powella oraz metodę gradientu sprzężonego Fletchera—Revessa. 


PRZYKŁAD ZASTOSOWAŃ 


Przedstawiony algorytm praktycznie zastosowano w procesie modyfikacji 
właściwości dynamicznych napędu głównego frezarki FWC-25 [9]. Problem 
modyfikacji właściwości dynamicznych produkowanej obrabiarki występuje 
w praktyce projektowej w przypadku korzystania z istniejących rozwiązań 
konstrukcyjnych przy projektowaniu nowej obrabiarki lub rodziny obrabia- 
rek. Obliczenia modyfikacyjne przeprowadzono zgodnie z algorytmem przed- 
stawionym na ryc. 2. Do oceny jakości dynamicznej rozważanego napędu 
przyjęto dwa wskaźniki [2]: 


U, = Ż e, (0,;) (31) 
oraz 
U;= | g (0)ez(o)do x) _g*(oj)e, (32) 
gdzie: 
Pw (Wy; A : e A 
e; (0) = Dirk dynamiczna podatność skrętna napędu mierzona na 
J 


wrzecionie względem momentu wymuszającego działającego na j-tą masę 
modelu napędu, 

e„(w) — podatność dynamiczna wrzeciona na kierunku p, 

g(wo) — funkcja masy, określona przez wartości współczynników rozktuti; 
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pewnej w szczególny sposób zdefiniowanej siły skrawania, uwzględniającej 
różnorodność stosowanych narzędzi i warunków obróbki. 

Na podstawie analizy rozwiązania konstrukcyjnego i kinematycznego 
oraz wyników badań analitycznych i doświadczalnych [10] ustalono, że dwie 
prędkości obrotowe N = 125 obr/min i N = 380 obr/min są reprezentatywne 
pod względem właściwości dynamicznych dla napędu frezarki FWC-25. Dla 
obu prędkości obrotowych opracowano dyskretne liniowe modele drgań 
skrętnych (ryc. 5) o postaci (4). Wyznaczone analitycznie strukturalne padra- 


Ryc. 5. Model fizyczny drgań skrętnych napędu głównego obrabiarki 


metry modeli poddano identyfikacji w oparciu o charakterystyki częstotli- 
wościowe drgań skrętnych napędu, mierzonych na. końcówce wrzeciona. 
Identyfikację przeprowadzono zgodnie z algorytmem przedstawionym w 
pracy. W wyniku tej identyfikacji otrzymane wartości strukturalnych para- 
metrów modeli zestawiono w tabeli 1. 

Dla modeli o zidentyfikowanych parametrach obliczono charakterystyki 
częstotliwościowe skrętnej podatności dynamicznej wrzeciona, na podstawie 
których wyznaczono wartości wskaźników oceny U, i U, dla obu prędkości 
obrotowych (dla N=125 obr/min U; jas = 3,241:10 %; U;js= 


Tabela 1. Parametry modeli drgań skrętnych napędu głównego frezarki FWC-25. 


J [da kg m”] k [da Nm] h [da Nms] 


13074: 5,0378 5,66573, 21,44 
PARA Z 2825130 4,26557, 8,91 
POM 332449 4,68983, 10,98 
94526 6,45067 2,985145 10,19 
29435 1,4788 2,341935 29,2 
4,4865 1,6259, S245 1:30 89,58 
LOST 1,8386 3,4965 1012,67 
6,52870' 2;560lEte 17482$5 141,16 
2,0497 7 7,1854$0 SACZAW EA 38,61 
24345 1,1592, 4,1852, 62,29 
200 2,09355 8,4248, 112,48 
C0. ke 7,1009ż6 314427 39,13 
3,925230 216,28 
7,901635 435,38 
2,9489, 67,11 


1 
2 
3 
4 
5 
6 
+ 
8 
9 
10 
11 
JE: 
13 
14 
15 


88 


= 1,158-1077; dla N = 380 obr/min U; zgo = 6,151:10"*; 


Us 2408005 10777. 

Z punktu widzenia przyjętych wskaźników oceny i ze względu na często- 
tliwości występujących wymuszeń roboczych najbardziej istotne są właści- 
wości dynamiczne napędu w pasmach częstotliwości <12, 200> Hz przy 
N =125 obr/min i <9, 180) Hz przy N = 380 obr/min. W pasmach tych 
obliczono wrażliwość charakterystyki amplitudowej przemieszczeń skrętnych 
wrzeciona względem współczynników bezwładności J i sztywności k. 

Analiza wartości współczynników wrażliwości względnej i logarytmicznej 
dla każdej z rozpatrywanych prędkości obrotowych pozwoliła określić para- 
metry strukturalne o największym wpływie na ocenę właściwości dynami- 
cznych. W przypadku prędkości obrotowej N = 125 obr/min są to współ- 
czynniki bezwładności o numerach 2, 5, 4, 1, 3 i współczynniki sztywności o 
numerach 2, 3, 4, 5, 1. Natomiast w przypadku prędkości N = 380 obr/min 
— współczynniki bezwładności o numerach 1, 12, 2 i współczynniki sztyw- 
ności o numerach 2, 12, 5, 3, 1. Wymienione parametry strukturalne są 
funkcjami parametrów strukturalnych, kinematycznych, materiałowych. 
Przeanalizowano zatem istniejące rozwiązanie konstrukcyjne z punktu widze- 
nia możliwości dokonania zmian przy założeniu minimalnych kosztów. W 
wyniku tej analizy ustalono, że zmianom podlegać będą cztery parametry: 


p+ — moment bezwładności koła zamachowego, 
pa — moment bezwładności koła pasowego, 

pa — sztywność przekładni pasowej, 

pą — sztywność skrętna wałka odboczki. 


Dopuszczalne zakresy zmienności tych parametrów wynikające z ograniczeń 
konstrukcyjnych podano w tabeli 2. 

W celu znalezienia rozwiązania optymalnego wykonano eksperyment 
numeryczny o liczbie doświadczeń K = 15. W każdym punkcie planu obli- 
czono wartości wskaźników U; i U, które następnie poddano normalizacji. 
Obliczono wartości kryterium globalnego Ię, które przyjęto w postaci: 


Tabela 2. Dopuszczalne zakresy zmienności parametrów będących zmiennymi decyzyjnymi 
w procesie modyfikacji 


Wartości parametru | Graniczne wartości parametru 
w istniejącym 
rozwiązaniu 


Wpływ na parametr 
strukturalny 


8,710 


1,246 


426,55 426,55 


4,36023, 4,36023, 4,4105%% i k; dla N= 125 
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1 4 
IG — 4 ). U, (33) 
i=l 


Analiza obliczonych wartości kryterium globalnego wykazała, iż optymalny 
jest wariant odpowiadający punktowi planu o numerze 11. Dla punktu tego 
wartość kryterium globalnego 1,,, =0,2 (największą wartość kryterium 
globalnego IG, = 0,69 uzyskano dla punktu planu o numerze 5). Wartości 
parametrów określających rozwiązanie optymalne są następujące: 


=. : RES 4 . od = 
pi, = 218-107%; pa, =8,54:107%; ps, = 610. 


W pracach [3—5] przedstawiony jest przykład kształtowania właściwości 
dynamicznych belki z własnym napędem frezarki FWD 25U, przy zastosowa- 
niu w algorytmie modyfikacji metody wariantów konstrukcyjnych. 


PODSUMOWANIE 


Przedstawiony algorytm kształtowania właściwości dynamicznych może 
być uważany za operator odwzorowania odwrotnego A ', gdyż umożliwia 
uzyskanie efektywnego rozwiązania zadania syntezy. 

Zakres zastosowań algorytmu nie jest ograniczony tylko do rozwiązywa- 
nia zadania kształtowania właściwości dynamicznych obrabiarek. Algorytm 
może być z powodzeniem stosowany przy rozwiązywaniu zadań syntezy 
bardzo szerokiej klasy maszyn 1 mechanizmów. 

Algorytm modyfikacji w połączeniu z algorytmem identyfikacji charakte- 
ryzuje się dużą dokładnością obliczeń. Uzyskane przy użyciu tego algorytmu 
optymalne wartości parametrów opisujących konstrukcję odznaczają się dużą 
wiarygodnością. Prawdopodobieństwo uzyskania tzw. „trafionej” konstrukcji 
jest bardzo duże. Liczba poprawek, jakie zwykle są niezbędne, aby opraco- 
wywane urządzenie spełniało założenia techniczno-ekonomiczne, ulegnie 
wówczas radykalnemu zmniejszeniu. Można dzięki temu uzyskać pokaźną 
obniżkę kosztów ponoszonych na projektowanie i wykonanie prototypów 
oraz znaczne skrócenie czasu niezbędnego na przygotowanie produkcji 
opracowywanej konstrukcji. 


STEFAN BERCZYŃSKI, KRZYSZTOF MARCHELEK 


SCIENTIFIC BASES FOR FORMING DYNAMIC PROPERTIES 
OF MACHINE TOOLS 


Summary 


In the paper the method of forming dynamic properties of machine tools is presented. The 
method is based on the interdependance of modelling issues, identification and modification. The 
identification of the models” parameters can be carried out on the basis of incomplete 
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information derived from experiments. Gradient methods are employed to minimalize the 
criterion of identification. The algorithm: of modification enables estimation of solutions accor- 
ding to different criteria. As decision variables in the process of modification those construction 
parametres are assumed which have the greatest influence on the accepted criteria of estimating 
dynamic properties. The optimal values of decision variables are calculated by the method of 
systematic scanning of points in the plan of the numerical experiment. The presented method of 
forming the experiment plan enables uniform distribution of the predetermined number of points 
within the domain of constraints. An example of utilization of the worked out algorithm of 
forming dynamic properties is presented. 


CT9DAH BOPYBIHbCKU, KIIHMLITOOĆ MAPXDJIEK 


HAYUHBIE OCHOBBbI OIIPEJNEJIEHMA CBOACTB 
ĄAMHAMMUECKMHX CTAHKOB 


Pe3oMe 


B paóoTe rpercTaBJIeH METOJ1 OINpEJHEJIeHHA CBOŃCTB JJUHaMMYECKHX CTAHKOB. OH orliepa- 
TCA Ha B3aMMO3ABHCHMOCTb BOIIDOCOB MOJIEJIMDOBAHHA, MIIEHTHDHKAIMM M MOJKQHKALINK. 
MneHrupukanuło IlapaMeTpoB MOJeJleń MO%HO IIDOBECTH HA OCHOBE HEIIOJJIHOŃ UH(DOpMaHMHH 
H3 OIIbITHbIX HKCCJIENOBaHHŃ. JIJIA MUHUMAJIM3ALHKH KPHTEeDHA KHEHTHQHKALNHM HPHMEHAIOTCA 
MeTOJ1bi TDANHEHTOB. AJITOpHTM MOJMbUuKa1MM MaćT BO3MOXHOCTb OLEHHTb peLLeHHE IIO 
MHOTHM KpHTepuaM. 3a IIepeMeHHble peLIeHHA B IIponecce Mojrmbukaliu OepyTca Te KOHCTPYK- 
THBHblie IIapaM€Tpbl, KOTODPbie OOJIbIIIe BCETO BJIMAIOT HA IIDUHATbIe KPHTEpHMH OLNCHKH JIMHAMH- 
ueCcKHX CBOŃCTB. ONTUMAJIbHbIe 3HA4€HHA PpELIAIOINHX IIEpeME€HHbIX OITIPEJEJIAŁOTCA METOJIOM 
CHCTEMATHUECKOTO IIEpECMOTPA TOYEK ILJIAHA UMCJIOBOTO JKCIIEDUMEHTA. IIpern1oxeHHbiń MeTO1 
(bopMMpoBaHHA NJIAHOB JKCIIEDUMEHTA JIAGT BO3MOXHOCTbŁ paBHOMEpHO pa3MECTHTb 3ANAHHOE 
KOJIMYECTBO TOUEK Ha 3ANaHHOŃ N.IOIHaNu. IlpencTaBJieH nupuMep npakTHYeCKOTO MCIOJIB30- 
BaHHA pa3paGOTaHHOTOo aJIropuTMa (b©opMMpoBaHHA JUHaMMHUECKHX CBOŃCTB. 
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